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ABSTRACT 


Identification and control of non-linear dynamical systems is a very complex task which 
requires new methods of approaching. This research addresses the problem of emulation and 
control via the use of distributed parallel processing, namely artificial neural networks. Four 
models for describing non-linear MIMO dynamical systems are presented. Based on these 
models a combined feedforward and recurrent neural networks are structured to emulate the 
dynamical system. Further, a procedure to emulate multiple systems in a single network is 
suggested. A method for finding a minimal realization of a network is introduced. The 
minimization greatly reduces the complexity of the network without degrading the operating 
performance of the network. 

This work also examines the application of artificial neural networks for adaptive control. 
The multiple-system approach is used to find an adaptive neural network controller for non- 
linear MIMO time-varying system in a direct model reference control scheme. The controller 
network is trained using a procedure called back-propagation through the plant, which was 
extended in this work. 

The application of neural networks is demonstrated on a longitudinal model of the F/A-18A 
fighter aircraft both with the undamaged aircraft and with a damage mechanism as a time- 


varying MIMO dynamical system. 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not have 
been exercised for all cases of interest. While every effort has been made, within the time 
available, to ensure that the programs are free of computational and logic errors, they cannot 
be considered validated. Any application of these programs without additional verification is 


at the risk of the user. 
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I. INTRODUCTION 


In many engineering and scientific applications a dynamical system or process of 
unknown structure acts upon measurable input excitations to produce observable output 
signals. One way of gaining knowledge of the system is by analyzing available data 
using a chosen candidate model to represent the system and applying an algorithm to 
adaptively change the parameters in the model, to cause it to function similarly to the 
unknown system when driven by the same inputs. The procedure is called system 
identification [LL87]. Controlling unknown systems is approached mainly by robust 
control and adaptive control, both relying on linearized or linear system theory to 
compensate for changes in the system. While robust control stays in the linear regime 
based on state feedback, adaptive control makes use of identification techniques to find 
a controller to compensate for parameter variations in the system as well as in the 
selected model, and regulates the output of the system according to some linear 
reference, resulting in a non-linear controller [AW89]. 

This work explores the subjects of identification and control of non-linear 
dynamical systems from aspects of artificial neural networks. This chapter describes the 
research objectives and contains a survey of previous work in the area of utilization of 
neural networks regarding dynamical systems. This chapter also summarizes the main 


contributions of the research and concludes with the organization of the dissertation. 


A. RESEARCH OBJECTIVE 

The research focuses on utilizing artificial neural networks with non-linear time- 
varying multi-input multi-output (MIMO) dynamical systems. The objective is to emulate 
multiple dynamical systems in a single neural network and to apply the procedure to form 
an adaptive controller for a changing plant. 

The F/A-18 fighter aircraft was chosen to represent a MIMO dynamical system. 
The undamaged and impaired aircraft represented multiple models, for which an adaptive 


model-reference neural controller was designed. 


B. A BRIEF SURVEY OF NEURAL NETWORKS INVOLVING DYNAMICAL 

SYSTEMS 

Recently, there has been tremendous interest in exploiting artificial neural networks 
(ANN) to perform tasks of system identification and control [MSW90], [SM92]. The 
identification is done by emulation of the dynamical system without having to find 
specific parameter values of a hopefully matching model. Techniques from linear and 
adaptive control were used to select a proper network structure for the task, as was done 
by Narendra and Parthasarathy for emulation of nonlinear systems [NP90], by Scott and 
Collins {SC90] and Bertrand and Collins [BC92] who used aircraft models for emulation 
and control of linear systems, and by Barto, Sutton and Anderson for control of nonlinear 
plant of inverted pendulum [BSA83]. Controlling took shape in different ways, usually 
with single-input systems as like Guez, Eilbert and Kam did for second order linear 


system [GEK88], or as Ha, Wei and Bessolo did for linearized second order 


approximation of the F16 aircraft. The works of Kraft and Campagna, who used a 
memory addressable type neural network controller [KC90], and of Nguyen and Widrow 
for the truck backer-upper using sequential emulation and control learning [NW90] in 
nonlinear systems, are other examples of single-input systems. 

For emulation and control the general approach has been to have a full state output, 
as Goldenthal and Farrell [GF90] and Anderson [AC89] did for inverted pendulum, as 
did Specht in his general regression neural network [SD91], and as was done in most of 
the works mentioned earlier. The approaches to control varied from neural network gain 
scheduling of classical controller of an aircraft by Burgin and Schnetzler [BS90], to 
finding an inverse plant such as was done by Bertrand [BD91] and by Levin, Gewirtzman 
and Inbar in a binary type inverse controller [LGI91], to self tuning or model reference 
adaptive control as was done in [KC90] and by Saerens, Soquet, Renders and Bersin1 
using back-propagation through the plant training [SSRB90]. 

The field of multiple multi-input multi-output systems requires further research and 
exploration. It is the purpose, subject and substance of this work to expand the 


knowledge and understanding of such systems. 


C. THE CONTRIBUTIONS OF THIS WORK 

The research defines an approach to emulating non-linear MIMO dynamical 
systems using artificial neural networks, as well as the limitations associated with the 
systems and networks. Furthermore, the concepts developed for emulation are used to 


find a model-follower neural network controller for a MIMO system. 


The research addresses the issue of treating multiple models in a single neural 
network associated with emulation and adaptive control of general dynamical systems 
from acommon class. A general scheme is presented for training a single network with 
multiple processes sharing a common dynamical structure and the same order of input 
and output. Additionally,a method for finding minimal realization of an artificial neural 


network, using singular value decomposition of the weight matrices, 1s developed. 


D. ORGANIZATION OF THE DISSERTATION 

Dynamical systems in association with neural networks are addressed in several 
chapters of this work. The disclosure of what a neural network is and how it is trained, 
is described in chapters II and II. The architecture, structure and behavior of artificial 
neural networks are presented in chapter II. The back-propagation algorithm used in this 
work to train the networks 1s presented with a few variations and extensions in chapter 
Iii. Chapter [IV characterizes the classes of dynamical systems with the associated 
structural aspects and training of neural networks used to emulate the systems. Chapter 
IV also present the approach to handling multiple models in a single network. A brief 
overview on adaptive control with implementation of control using neural networks is 
given in chapter V, which also establishes the procedure of reconfigurable neural-net 
controller. After having an understanding of what artificial neural networks are, how 
they are trained and how they specifically apply to the problems of emulation and control 
of MIMO non-linear dynamical systems, some design considerations for dealing with 


dynamical systems are highlighted in chapter VI. Chapter VII describes the F/A-18 


aircraft as the dynamical system chosen to demonstrate the procedures which were 
developed in the previous chapters and is followed by a presentation and a discussion on 
the results in chapter VIII. Chapter [X closes with conclusions and recommendations for 
further research. To summarize, the research shows how identification and adaptive 


control problems can be solved quite straightforward using neural networks. 


HW. NEURAL NETWORK THEORY 


A. INTRODUCTION TO NEURAL NETWORKS 

A neural network, or better called artificial neural network, is a parallel, 
distributed information processing structure consisting of processing elements that are 
interconnected among themselves via unidirectional connections. The artificial neural 
network (ANN) has the ability to learn from experience, to perform abstractions of inputs 
with relevant information, and to generalize its knowledge from previous results. The 
massively connected network is highly parallel because the operations of each processing 
element occur simultaneously. Associated with each connection is a weight. The 
connection weights are the dynamic elements that change when the network is trained. 
The goal of the training is to adapt the weights so that the network performs in a desired 
fashion in response to a given input. 

The architecture of artificial neural-network has three main levels: topological, data 
flow, and neurodynamics. The architectural and algonthmic factors that compose any 
neural network are: the interface with the environment, the propagation paradigm, the 
learning rule and the recall or test phase. The structure and behavior of artificial neural 
networks are set by the layers that are built of processing elements and the activation 
function of each processing element. This chapter describes the architectural and 


structural factors and, where appropriate, gives a mathematical description. Three basic 


categories of networks that find use in this work are then described and are then followed 


by a discussion on the existence of networks suitable for this research. 


B. ARTIFICIAL NEURAL NETWORK ARCHITECTURE 

The architecture of ANN is described in three levels: topological, data flow, and 
neurodynamics. The presentation here will follow the guidelines of Neural Computing 
by NeuralWare, Inc. [NC91], who developed the basic software for this research. 

The topological level defines the structure of the network and its interface with the 
environment. The topology also describes how processing elements (PEs) are grouped 
into layers and how PEs are interconnected to one another. A /ayer here is defined as 
a group of PEs that act in a similar fashion. In Figure 2.1 a basic schematic of an 
artificial neural network is presented, consisting of an input layer, an output layer, and 
one hidden layer. The circles denote the processing elements and the lines denote 
connections. The term hidden layer refers to a layer which does not have a direct 
connection with the environment. The hidden layers are the essence of the network and 


extract the features that enable the network to solve nonlinear problems. 
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Figure 2.1: A Basic Schematic of an Artificial Neural Network 


The second level involves data flow and describes how information is presented to 
the network, is passed from one processing element to another, is acted upon within a 
processing element, and is delivered to the outside world. The data flow level defines 
the learning rule and the recall or testing process. 

The third level, neurodynamics, describes in detail the operations that act upon the 
data within a processing element. This level defines the functions and the behavior of 
the processing element. 

The following paragraphs describe the architectural factors in more detail. The 
description will start with the environment in which the network operates and will end 


with a scrutiny of the inner structure of a single processing element. 


1. Environment 

Recognizing and understanding the environment in which the neural network 
has to perform is as important as choosing the network itself. Most of the time the 
environment is the factor that dictates the choice of the neural net, its architecture, and 
its learning scheme. It is tmportant to have a proper interface, including scaling and 
normalizing the data. In particular, the I/O, its timing and characteristics, determines 
the number of inputs, outputs, and other criteria such as convergence and 
pre/post-processing. 

After the proper interface on the network is selected, the inner-structure of 
the network that propagates the data from the input layer to the output has to be chosen. 


The structure is governed by the propagation scheme. 


2. Propagation Scheme 
The general structure of the network is controlled by the propagation scheme. 
The propagation scheme defines how processing elements are connected to one another. 
Figure 2.1 is an example of connections in a fully connected feed-forward network, from 
the input elements, through the hidden layer, to the output nodes. The propagation 
scheme governs the connectivity matrix structurally and by type. The connections may 
be sequential, by layer, random in order and occurrence, or they may be dictated by the 
properties of the specific element. 
Now that the network is formed, it has to be trained to perform the desired 


task. The training is done via a learning rule. 


3. Learning Rule 

Learning is the phase in which the connection weights are adapted according 
to relation in the input data to some other process. The kind of process categorizes the 
learning type into one of three categories: supervised learning, unsupervised learning, 
and reinforcement. 

Supervised learning is a process where a desired output corresponds to a 
given input. Both, input and desired output, should be supplied in the learning phase, 
during which the weights are updated via some learning rule to match the input-output 
pair. When the input-output pair is the same vector, the network is referred to as being 
auto-associative. If the desired output and the input vectors are unequal the network is 
called hetero-associative. The work of this este uses supervised hetero-associative 


learning. 


Unsupervised learning occurs when no desired output is presented and the 
network groups the input presentation into some meaningful clusters. 

Reinforcement is a process that falls between the two above-mentioned 
categories. In this case, there is an indicator, a teacher, that indicates whether the 
response to the presented input is satisfactory. 

In any learning category the process is governed by a certain learning rule. 
This learning rule dictates the change of connections’ values or weights to achieve the 
goal of the training. The parameters that control the learning rule may change while 
training. The learning phase may need thousands or even millions of presentations until 
convergence to a satisfactory result is achieved. The learning sequence usually involves 
a feedback process in which the weights are updated. The specific learning rule that was 
used in this research, called the ee algorithm, is described in detail in 
Chapter M1. 

When the network is trained to satisfactory degree the learning is stopped and 
the weights in the network are frozen. The network can then be used to perform the task 
it was trained to. This implementation is governed by the recall phase which is described 
next. 

4. Recall Phase 

The recall or testing phase is represented as a feed-forward pass in the 
network that gives a response to a certain presentation at the input layer. The recall pass 
might be a subprocess of the learning phase in which the network output is calculated 


before the adaptation action is taken due to output error with respect to the desired 


10 


network response. In essence, after learning is completed, testing 1s the practical and 
useful phase. 

Recall phase is strongly related to the propagation scheme. To better 
understand the propagation through the structure, a basic multi-layer feed-forward 


network is discussed and analyzed next. 


C. ARTIFICIAL NEURAL NETWORK STRUCTURE 

The behavior of an artificial neural network (ANN) is determined by the network’s 
global and internal structure. This section gives the mathematical description of behavior 
of an ANN. 

The nomenclature adopted here and throughout the dissertation is as follows: 
Scalars are denoted in upper-case letters. Vectors are denoted in bolded lower-case 
letters. Matrices are denoted in bolded upper-case letters. Elements of vectors or 
matrices are denoted with corresponding unbolded lower-case letters with appropriate 


indexing. Additional indexing, when required, is given and explained in the text. 


1. Multi-Layer Neural Network Structure 
Multi-layer networks belong to a family of networks consisting of an input 
layer, output layer and one or more hidden layers. Figure 2.2 shows a basic three layer 
feed-forward artificial neural network and its block diagram representation. The 
presented network is used for the analysis of the structure and operation. A basic 
network was selected to simplify description and understanding without the loss of 


generality in the more complex multi-layer networks. The network shown is fully 


LI 


connected between consecutive layers and a bias term is connected to every processing 


element. 
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Figure 2.2: A Three-Layer Feed-Forward Artificial Neural Network 


The network is built from layers that take the same general form. 


Mathematically a layer is described as 
y = Mx) eee 


where JN, is a function that maps x€R"° to yE&R™. For a specific layer, r, the 


nonlinear function N,€ R°—R™ is written as 
N(x) = F[W"x +5" (2:2) 
where F, is a diagonal nonlinear operator which represents a set of m nonlinear 


activation functions f€R-R that are identical for the current [r]" layer [NP91]. The 


operand of F is the weight matrix W’'€ R"*™ times the input vector x © R" summed over 


12 


all connecting nodes and with a bias term BER". A multi-layer network containing R 
layers, excluding the input layer, is denoted as W,,,,,__,,- Each [{r]” layer consists of 1, 
nodes (processing elements). The network, N, transforms an input vector x€ R™ to an 
output vector y€ R™ propagating via (R-1) hidden layers, each consisting of 71,,7,,...,"g, 
nodes, respectively. 

For the three-layer feed-forward network shown in Figure 2.2, the output of 


the network, y, as a function of a given input vector, x, is given by 
y¥ = Ninaa(X)=F,{Welr,[WOF (Wilx +b!) +b21 +58} = NN,N,(x) (2-3) 


The input vector, x, is presented to the network and processed in parallel in each layer. 
The information ts passed in cascade from layer to layer and presented at the output layer 
as the network’s outcome. 

Each layer in the network is built from processing elements that work in a 
similar fashion. The processing elements carry the processing task of each layer. 


Another step down to the basic element of the ANN is taken next. 


2. Processing Element 
The processing element (PE) is the basic element in any artificial neural 
network. It has virtually as many inputs as desired and one output. The inputs are 
summed to create net input for the given PE. A bias is added and the result is acted 
upon via a scalar activation function. The output of the activation function is the output 
of the PE. A diagrammatic representation of a processing element is given in Figure 


2.3. The basic PE configuration just described is used in the present research. 
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Figure 2.3: A Basic Configuration of a Processing Element 


The mathematical description of a processing element in the form presented 


in Figure 2.3 is 
yy ea) (2.4) 


[il = 3 Wi, sig (255) 
where y/” is the output of the PE mite [r]® layer and ” is the activation function 
that acts upon the summation value J/”. The summation is done over all the outputs .x, 
of PEs that are connected to the / PE in the [r]” layer via a weight w,/7 and the bias 
b{". It is noteworthy that throughout the work the first subscript of a weight component 
relates to the target PE whereas the second subscript relates to the source PE. 

Variants of a PE’s structure may include different summation functions, such 
as cumulative sum, maximum, minimum, product, etc. The output can either be passed 


directly or a competitive process between PEs of a given layer can be established 
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(Kohonen layer). The output value can feed any number of other PEs’ inputs, including 
its own. 

An important part in determining the output of a PE is the activation function. 
Selecting a proper activation function is crucial to successful network training. A more 


elaborate discussion Is given next. 


3. Activation Function of a Processing Element 

The activation function is a scalar function that acts upon the net sum of the 
effective inputs. In general, activation functions can be linear, semi-linear or nonlinear 
functions, and may be non-differentiable. A nonlinear network is achieved by 
introducing PEs with nonlinear activation functions. In linear processes a linear 
activation function can yield better results, but generally the more interesting cases are 
nonlinear. In some cases the weights can even be pre-predicted as will be shown in 
section C of Chapter IV. 

The two activation functions used in this work are the linear function, and the 
hyperbolic tangent. 


The linear function, given by 
f(D) = I+ Gain+b (2.6) 


Where / is the summation value multiplied by Gain parameter and 0 is a bias term. 
The hyperbolic tangent, a continuous monotonic squashing function, which 


maps the input into a range from minus one to one is given by 


I = I- Gain 
: ar ape (2.7) 
(stan) = ae 
rr, 
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The hyperbolic tangent is differentiable, and moreover, the derivative can be expressed 


in term of the function itself as 


FD) = (+f) * (1 -fO) —_ 
The differentiability condition is needed in the application of the back-propagation 


learning algorithm and an analytical differential as in equation (2.8) saves computations. 

Another function commonly used in networks trained via the back-propagation 
algorithm, and which is similar to hyperbolic tangent, with the difference that it maps 
to a range from zero to one, is the sigmoid function. The sigmoidal function also shares 
the property that its derivative can be expressed in terms of the function itself. Figure 
2.4, below, shows graphically the functions and their derivatives. Some other activation 
functions that are used in various neural networks include, but not limited to, sine, 


signum, step-function and piecewise linear. 
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Figure 2.4: Hyperbolic Tangent and Sigmoid Activation Functions and Derivatives 


Having finished the general structure of an ANN, the specific configuration 


of the multi-layer networks used in this research are descnbed next. 


D. CATEGORIES OF MULTI-LAYER NETWORKS 

Multi-layer networks are divided into two main categories. The first and most 
common is the feed-forward network, as shown earlier in Figure 2.2. The second is the 
recurrent network. Introducing the recurrent network is essential when dealing with 
dynamic systems. A block diagram of a recurrent network is shown in Figure 2.5. A 
multi-layer network may include either or both types. A block diagram of a basic unit 
of a generalized multi-layer network is shown in Figure 2.6. The principle of combined 
recurrent and feed-forward network can be further extended to include information from 
prior input and output vectors, as shown in Figure 2.7. The idea of using general 
networks in conjunction with dynamic systems, and the detailed structure of the block 


denoted by N, will become clear in Chapter IV. 





Figure 2.5: A Recurrent Figure 2.6: General Multi-Layer Building Block 
Network 
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Figure 2.7: General Regression Network 
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E. ARTIFICIAL NEURAL NETWORKS AS GENERAL APPROXIMATORS 

Before choosing a network with a specific structure one may ask if there is a 
network that can do the desired task. The answer to such a question Is in the affirmative 
and is based on two major theorems. The theorems express the existence of artificial 
neural networks as function mapping devices. ANN is found to be capable of 
approximating functions from a finite dimensional space to another, to any desired degree 
of accuracy with only one hidden layer. 

The first proof of the capability of a feed-forward multi-layer network to represent 
continuous functions of many variables by superposition of continuous functions of one 
variable and addition was given by Kolmogorov’s superposition theorem. Hecht-Nielsen 
entitles the theorem as "Kolmogorov’s Mapping Neural Network Existence Theorem" 
[HN9I1, pp. 122-124]. The theorem shows the existence of a neural network with linear 
input and output layers, with as little as one hidden layer, as being capable of mapping 
any continuous multi-vaniable function. A corollary gives an upper limit to the number 
of hidden nodes needed for the mapping (27)+1, where Mp is the number of input nodes). 

However, for every mapping task a different unknown activation function is needed. 
Other expansions of this theorem show the existence of neural networks as mapping 
devices on different but very wide and practical function spaces [HN91, pp. 131-133]. 

An important theorem based on Stone-Weierstrass’ theorem, given by Homik, 

Stinchcombe and White, shows that ANN are general approximators of a bounded 


function using an arbitrary squashing activation function, i.e., a sigmoid or a hyperbolic 
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tangent in an unknown but large enough number of PEs to get the desired accuracy 
[HSW89]. The results show that multi-layer neural networks are universal approximators 
in mapping any measurable function to a desired accuracy and any lack of success in 
application must arise from inadequate learning, an insufficient number of hidden PEs, 
an inappropriate activation function, or the lack of a deterministic meaningful relationship 
between the input and the desired output. 

However, neither theorem really helps to build a “correct” structure of the 
network, although both show that three layers are enough. Since three-layer network 
approximators may require an impractical large number of hidden processing elements, 
in real world problems multi-layer (>3) networks are used to find a tractable solution. 
The theorems, although not helping in the selection of a proper network structure, do 
give confidence that an appropriate network exists. Thus, although a suitable network 
exists, there 1s no guarantee on the existence of a learning rule to train the weights. 

An understanding of the physical nature of the problem may assist in choosing a 
“good” number of hidden layers and hidden nodes. Until the present, choosing network 
parameters was a kind of art, usually done by tral and error. A more mgorous way of 
selecting the number of hidden layers and how many nodes each layer should include is 
presented in this dissertation as minimal realization of a network (Chapter VI.D). 
Although the method cannot predict in advance a good choice of number of layers and 
PEs, the rationale for going through the procedure of finding minimal realization of the 
network is to achieve the same accuracy with less network complexity. Reassured that 


a suitable network exists, a suitable training algorithm needs to be chosen. The 
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algorithm used in this work, which has been successfully applied to numerous multi-layer 


feed-forward networks, is the back-propagation algorithm, described in the next chapter. 
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I. ANALYSIS OF BACK-PROPAGATION ALGORITHM 


Among the many types of learning algorithms, the back-propagation algonthm has 
been chosen for this research. This learning rule is also called the generalized delta rule. 
The name back-propagation is based on the way the errors are handled in changing the 
weights during learning and is attributed to Rumelhart [RM86, vol. 1, pp. 319-364]. 
This algorithm is implemented on a hetero-associative network. The standard back- 
propagation learning rule is derived, followed by four common variations. An extension 
of the standard algorithm is also presented, with some new highlights used in this 


research. 


A. ARCHITECTURE 

Multi-layer networks based on the back-propagation algorithm are feed-forward 
networks. They consist of an input layer, an output layer, and at least one hidden layer. 
A typical fully connected feed-forward network with two hidden layers is presented in 
Figure 2.2. To keep the discussion general, the rest of the derivation will not refer 
specifically to a particular number of hidden layers or hidden nodes. The representative 
network shown in Figure 2.2 is merely for demonstration and ease of explanation and 
does not affect the generality of the back-propagation algorithm. 

The back-propagation learning algorithm is used to train multi-layer networks to 
map or approximate a process or function y*=.¥ (x); # ER™*", represented by input- 


Output pairs Peavale, where d denotes desired values and p stands for different 
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presentations. An analytical representation of the function may not be available. After 
the training is complete to a desired accuracy the weights are fixed and the network then 
generates an approximation of the desired function, y, fora given input, x, as y=A(x). 
The rationale for this method of adaptation is given in the next section and the basic form 
of the learning rule is denved. Many variations have been suggested, mainly to increase 


convergence speed and to avoid reaching local minima. 


B. MINIMIZING GLOBAL MEASURE 

The goal of the learning process is to minimize a global measure by modifying the 
weights in the network. This measure is usually taken as a function of the error at the 
output of the system. The global measure can in some cases be interpreted in terms of 
the energy content of the system. Minimizing a global quadratic error is well treated in 
the literature and is based on the gradient descent method, sometimes referred as the MIT 
rule, e.g., [BH75, pp. 237-240] or [LL87, pp. 282-288]. The discussion here will 
concentrate on the discrete domain where a continuous domain equivalence also exists 
[NP91]. 

Let a measure, J, as a function of the parameter space, 6, be defined as 

N 
J(9) = ie L[e(0,k)] (3.1) 

where & is the discrete time indicator and N is the number of samples taken. Among 
many sensible error functions, L, that can be used, a commonly used one that acts on 


a network’s output error, e, is the global quadratic error, defined as 


Lle(6,k)] = ~||e@,4)||? = 5e7e (3.2) 
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Minimization of the measure J with respect to the parameters @, is done by going 
in the negative direction of the gradient of J in the parameter space. Hence, the change 
in the parameters is proportional to -V, J such that, 

6,= 6,+48 , and (3.3) 


A8 = => Vi (3.4) 


where the gradient is defined as V,J4 ie , and 7 1s a constant of proportionality. 
Utilizing the gradient descent method to minimize J, the gradient is determined from 


equation (3.1) using the chain rule 





_ ry Ble (0,4) dL[e(9,k)] _ de(9,k) (325) 
<~ —9e(0,k) 06 


For optimizing dynamical systems best results are achieved when Noo, but in practice, 
N 1s chosen over a finite interval, 7, [k-7+J1,k]. In static processes N=1 is sufficient 
[HN91, pp. 136] and usually is the case implied in the back-propagation algorithm. 
Using the global error function as in equation (3.2), letting i represent the instantaneous 


time, the performance measure takes the form 


k 





J(0,k) = = SY) e(0,i)7e(0,i) , (3.6) 
isk-T+l 
and the gradient is 
k é 
VJ =13 e(0,i)7 | 2G) (3.7) 
T ekoT+l 06 


The parameters, 9, are then updated every 7 samples. A single parameter, 6,, at time 


t=(k+1),is updated as 


fi 


e(0.i) (3.8) 





k 
- _ Be a _ I de(0,1) 
O(k+1) = O(k-T+1)  . 
i=k-T+ ] 
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The above equation implies that the gradient method can be implemented if the value 


9e(6,i) 
a6 


The major disadvantage of the gradient descent method is the slow convergence 


can be computed on-line for all values of i in the interval [k-7+J,k]. 


rate near the optimum [BH75, p. 238]. 
The next section shows how this kind of parameter optimization is applied to the 


back-propagation algorithm. 


C. BACK-PROPAGATION LEARNING RULE DERIVATION 


Suppose the network in consideration is from the class NM ney aS defined in 


no,nt... 
section II.B.5, where the input vector is x=[X,,x,,...,X,,], and the output vector is 
y=[),,¥2,--->Yng]. The parameters in the network are the weights W, which consist of 
all relevant w,/7; r=1,2,...,R, as defined in section II.B.6. The training sets are [x,y], 
for every presentation, p, and y’ has the same dimension as y in order to train all the 
outputs of the network. The network propagates the input vector to produce an output 
vector, y, which is compared with the desired vector, y’ and result in an error vector, 
e, where e,=y,,7- y,; m=1,2,...,N,g. If there is no difference, no learning takes place. 
Otherwise, the weights are changed in the direction of the negative gradient, -V,,F, to 
reduce this difference, i.e., the method of steepest descent. The performance measure 
is defined here as E, in correspondence with the nomenclature used in the neural- 


networks related literature, corresponds to the total-error energy in the system, in the 


mean square sense. In accordance to equation (3.6) the measure is given as 
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N 
al 8 
a er (3.9) 


je 5 ele = > (ya-y,)? (3.10) 
where E,, is the error energy ina single presentation p, over all the n, output PEs. The 
dependance of the measure, EF, on the weights, W, is through the output y. The change 


in the weights as in equation (3.4), is proportional to the negative gradient, 


_ OE 
aw 


_ PP yap sal JE (W) (3.11) 


[r} _ 
Aw; n a awit 


y 
p=l 








The coefficient of proportionality, 7, 1s called the learning rate. In the batch update, 
the number of presentation, N, 1s to be chosen according to the problem, and does not 
necessarily have to include all presentations. Furthermore, in practice a small finite 
number is taken, usually N=1, 1.e., an immediate update. 

For simplicity, without affecting the generality of the mathematical calculations, 
cascading fully connected layers will be assumed. For clarity, the governing equations 


for data propagation, equations (2.3) and (2.4), are repeated here. 
i" = pp wy (3.12) 
aft = pina) G.13) 
The value z, represents output of hidden PEs. At the input layer //”=.x,, and at the 


[R] 


output layer z;""=y,. The error at the output layer is given by 


e =y-y, (3.14) 
The gradient in equation (3.11) is derived assuming that N=1. The end result will 


be generalized for any N. 
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Since any functional dependence of E, on w,/” must be through //”, one can write 


the gradient using the chain rule as 
) a) Hl 





= ame (3.15) 
aw? ari” awh? 
Define a local error as, 
dE 
BE = f (3.16) 
Ge 


Substituting equation (3.16) and using equations (3.12) and (3.14) in equation (3.15) 


gives the negative gradient as 


r=-| 


OE = sll, r] Swe . ft. gill (3.17) 


Pp 
{r] 
ij 








[fr] iq “q 


dw, 


Ow:. g=1 





where z/”” is the output signal of the /” PE in the [r-1]" layer in the forward pass of the 
p” training presentation. 


Evaluating the local error term 6, and using equation (3.13) in the derivation gives, 





sil dE, _ OE, zi” _ OE, rrp ftly (3.18) 


where /f’ is the derivative of the activation function with respect to its argument, 


“ . eae dE 
evaluated at J/”. This illustrates the need for the existence of the derivative term g 





az! 
At the output layer z/*/= y,, then 
JE, OE, . a 
at ype” Hamas <_ 
git =¢ fer) (3.20) 
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At layers other than the output layer, i.e., hidden layers, use of the 


multidimensional chain rule gives 


ee [rel] 
aE, = dE, ali 3.21) 


af? get oh! agi” 
Using equations (3.12) and (3.18) in equation (3.21) yields 


* “Ya a. wht) and (3.22) 


a7 = pirat). Sal wie sical 

q=l 
In essence, no particular restriction was placed on the error function, E. The only 
requirement is that the activation function will be differentiable at least once. A 
recursive mechanism to back-propagate the error at the output layer, equation (3.20), to 


the hidden layers, equation (3.23), is formed. Altogether, it gives a straightforward 


learning rule, where according to equations (3.3), (3.8) and (3.11), 


We = Wy + Zw, (0) = Wj,” ¥ 0 (3.24) 
where, 
5 = fra: 3 ow pa vO0<r <R, (3.25) 
g=l 
and for the output layer, 
Se) = fh) - yy’ = y)) a, 
In the general case where N21, equation (3.24) becomes 
we = wy! + DD) eee C27) 


and equations (3.25), (3.26) are with addition of subscript p on all variables. 


IT 


Figure 3.1 [NP91, p. 254] gives a diagrammatic representation of the process 


including the forward and backward passes. 


a an. 
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Figure 3.1: A Diagrammatic Representation of Back-Propagation Architecture 


As shown, in the backward pass the first derivative of the activation function is 
used as well as the value of the function itself. For the sigmoid and hyperbolic tangent 
activation functions, in addition to their desired behavior emphasized in section II.B.7, 
their derivatives can be expressed in terms of the function itself. Representing the 
derivative by the value of the function itself saves in the number of calculations needed, 
and it may also enable a much simpler hardware implementation of a multi-layer, 


feed-forward neural network trained via back-propagation algorithm. For example, using 
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hyperbolic tangent activation function, equations (2.7) and (2.8), in equations (3.25), 


(3.26) gives 


al = (LedlNaaly Salt. wit! | and (3.28) 


q=l 
6" = (1+y)(1-y)- OY - yD (3.29) 


respectively. 


D. COMMON VARIATIONS ON THE STANDARD BACK-PROPAGATION 

ALGORITHM 

A difficulty with gradient descent-based algorithms, including back propagation, 
is choosing a suitable learning rate - 7. Updating the weights by means of equation 
(3.24) is based on the assumption that the error surface is locally linear with respect to 
the step size. At points in which the error surface is “highly” curved, too big a step size 
may cause the process to become unstable and diverge. On the other hand, keeping the 
learning rate small yields a very slow convergence, and the process may get trapped in 
a local minima. Many methods have been developed to overcome this phenomenon and 
create a more stable update process and much faster convergence. Some of the simple 


and effective methods are briefly summarized here. [NC91, pp. 94-99] 


NS, 


1. Momentum Term 
A momentum term is introduced in equation (3.27) which gives significance 


to the previous change in the weight. The equation with the momentum term is, 


wi” (k) = wi” (k-1) + Aw,” (k) 
(3.30) 


Awl(k) = 16k) (k) + wAwy(k-1) 


where 0<a<l is a scalar momentum term between zero and one. Introducing the 
momentum term enables faster learning with low learning rate, 7, and also allows an 


increase in 7 and @ before the learning process goes unstable. 


2. Cumulative Update of Weights 
Implementing equation (3.27), without the normalizing term ('/y), (with or 
without the momentum term) also improves the convergence rate, by the fact that A 
weights in the wrong direction are averaged out. 
This approach may lead to inadequate learning in some data sets. This is 
because for well-behaved learning, the training sets have to be well randomized. If the 
data comes from more than one source, i.e., multiple-function mapping, desirable 


changes might cancel each other. 


3. Fast Back-Propagation 
Another variation on the standard back propagation that greatly improves the 
convergence rate, although it is more sensitive to divergence, is the addition of the error 


to the PE value. 
Awit = -57-(2!7 + Bar") Co 
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where 0<6 <2 1s a factor, which when set to zero yields the standard back-propagation 
algorithm. It can be seen that asymptotically when 6-0 (as training progresses), this 
learning rule approaches the basic back-propagation algorithm. The momentum term as 


discussed in section 7 may also be applied. 


4. Modifying the Derivative 
As seen in Figure 2.5, at the extremes where the activation function becomes 
saturated, the derivative of the activation function goes asymptotically to zero, and thus 
the PE stops learning, since Aw«/’. By adding a small offset to #/, the PE continues 
learning even when it is saturated. This can in turn speed up the learning. 

Another significant variant of the back-propagation algorithm, called back- 
propagation through the plant, that treats the problem where the desired values at the 
output layer are not known directly, is presented in the next section. It should be 
mentioned that any of the four variants described above can be further applied to the new 


algorithm. 


E. BACK PROPAGATION THROUGH THE PLANT 

The following training algorithm was developed in the context of controlling a plant 
with only qualitative knowledge on the plant behavior by Saerens and Soquet [SS89]. 
An extension of the algorithm is proposed here to give a more accurate implementation 
of the gradient descent, and results in faster convergence. The change comes in 
estimating the Jacobian of the plant, and keeping the process within the domain, an 


important issue when dealing with non-linear systems. 
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1. Architecture and Formulation 

As was shown, in order to train a network, the pair of input and 
corresponding desired output values need to be supplied. There are instances where the 
neural network is cascaded to another physical process, and the ANN output is connected 
to the process or plant input, as shown in Figure 3.2. Suppose the network in 
consideration is from the class N,,,,,. ne: aS defined in section II.B.5, where the input 
vector is x=[X,,X,,...,X,,], and the network output vector is u=[U,,l,...,U,,]. The 
parameters in the network are the weights W, which consist of all relevant wi” 
r=1,2,...,R, aS defined in section I.B.6. The output of the network u, is used to excite 
a non-linear dynamic process, &. The output of the non-linear process is given by 
y=[),,y2---»Ys], Where S is the number of outputs of &. The training sets are [x, ¥1, for 
every presentation, p. The formulation of the combined system, with special attention 
to the indexing used in the derivation of the learning rule, is given by 


y(k+l) = ® [u(k), v(K)] 
u(ky = Nn ie (k), WCK)] 


3:32) 


where v(k) are internal parameters of @ that do not depend on the current input u(k), but 


may depend on previous inputs, W represents the weights in the network, and 


u(k) = {u,(k), l=1,...,n,} 


y(k) = ty (4), m=l,...,S} 


(3.33) 





Figure 3.2: Artificial Neural Network preceding a Physical Process 


a2 


The desired values are known at the output y, but not at the intermediate 
point u, which are the values needed for training a hetero-associative neural network. 
In cases like this, a variant of the back-propagation algorithm is used, called back 
propagation through the plant (BPTP). The idea is to transfer the error information at the 
output of the plant, through the plant, to the output of the network u. There are two 
ways of doing this. 

In the first a neural network is used to emulate the plant or process, and the 
weights are frozen after training. A new neural network is built and connected in 
precedence to the emulated process. The combined network is trained using regular back 
propagation. The error at the output y is transferred backwards through the fixed 
weights of the emulated plant. The transferred error is further back propagated through 
the preceding network. The weights in the preceding network are updated until a small 
enough error is achieved at the output of the emulated plant. 

The second method avoids the need to have a neural network representation 
of the process/plant. The error at the output of the plant y, is transformed to a 
corresponding error at the output of the network u. From that point it is possible to use 
the standard back-propagation algorithm. The procedure of back-propagating the error 


through the plant is derived in the next section. 


2. Derivation of Learning Rule for Back-Propagation Through Plant 
As before, the objective of the training is to minimize an error function that 


acts upon the output values due to changes in the weights. Define a measure on the 
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q 
output as in equation (3.10) as E, = 3 (iy, -y )’. The parameter space is W, where 
* m= 
w,/7 is the weight connecting the output of the /” PE in source layer [r-1] to the input of 
the i“ PE in target layer [r]. The development will follow the procedure set in 


subsection C. The change in the weights is proportional to the negative gradient of the 








measure with respect to the weights as in equation (3.11) Aw," oc Whol = -y at, 
Following the exact derivation as in equations (3.15)-(3.18), gives 
; ~s = 672°", and (3.34) 
ij 
gf = Eo. preg (3.35) 





az” 


The partial derivative of the measure with respect to the output of each PE needs to be 
evaluated. Starting at the output layer, using a multidimensional chain rule and utilizing 


equations (3.10), (3.11) and (3.32), noting that at the output layer z”!=u, gives 


aE, S. dE, a s 
_) eae (3.36) 
du, m=] dy. du, m=] ie 


where J, is the Jacobian of the function @ with dimensions Sx7,, and the (m,/) element 


is given by ae Substituting equation (3.36) into (3.35) yields 


Ou, 


s 
6 = parryi) SS e. Te (3.37) 
m=] 


and therefore, 
AW, (= 1 Cuca (3.38) 


which has the same form as equation (3.24) in the generalized data rule, and equation 


(3.26) is replaced by equation (3.37) as the back-propagated error at the output layer. 
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Following on with the derivation, equations (3.21)-(3.23) give the same results, 1.e., 
equation (3.25). 

The BPTP learning rule is a recursive procedure identical to the standard 
back-propagation algorithm, with the modification on how the error at the output layer 
is calculated. The same result under different derivation can be found in [SS89]. 
Recognizing in the classical back-propagation scheme that @(u)=[,-u, and J;=[,, 
where L,, is the identity matrix of size ng, equation (3.37) reduces to (3.26) exactly. 
This shows that the BPTP variant is a direct extension of the basic rule, or rather, the 
basic rule is a particular case of the former. Hence all the other variants described in 
subsection D apply to back propagation through the plant as well. 

In the derivation, the only parameter that does not have a closed from for the 
general non-linear case, is the Jacobian. Researchers treat it in different ways, and 


another approach to estimating the Jacobian is proposed next. 


3. Estimating the Jacobian of the Plant 

Close examination of equation (3.37) leaves the Jacobian J, a process 
dependent parameter, that not always has a closed form formulation. Jordan [JI89] 
suggested a way to calculate an estimation of the Jacobian using a neural network 
emulation of (uw). Saerens and Soquet [SS89] showed that if no details are available on 
the Jacobian, knowing its sign might be sufficient, given this is a pron known. Because 
this 1s a rough approximation, convergence is rather slow, since the descent is not in the 
direction of the gradient. In nonlinear systems, monotonic behavior is not guaranteed. 


In off-line procedures, where the process P(u) is simulated, if the exact value were 


oa) 


unknown or were difficult to calculate, a simple numerical differentiation may be used. 
A similar technique was used by Psaltis, Sidris and Yamamura [PSY88]. As was 
mentioned earlier, the exact value is not necessary [SS89], but a close value is preferable 
for more rapid convergence. Thus, an “influence coefficient” may be considered instead 
of the Jacobian. Assuming ® is well behaved at the point of calculation for a small 
enough Au, backward or forward differences are sufficiently accurate and simplify the 
calculations in computing an estimate of the Jacobian. The choice of backward or 
forward differences is based on the previous state of u in order to ensure staying within 
the domain. The influence coefficients are calculated at every iteration for every state 
in the vector u. Denoting the iteration sequence by the index &, the choice of difference 
method is as simple as 

if u(k) > u(k+1), use forward differences. 

if u,(k) <u(k+1), use backward differences. 
The calculation is done for each state in u uci leaving the rest constant. 

The influence coefficients are calculated as follows: 

Using forward differences 


D (Uj sly, +++ Ue tA Uy... Uy K)-P (Uj las «++ sly +» Uy »K) 


J _(k = is (3.39) 
mi) Au, 
Using backward differences 
h 53 ose (LER Ub SIC) DNC nn sd tet, 
7b _ (ltt u, U, k)-® (u,,u, as u, k) (3.40) 


AU, 
Where ®,(- ,k)=y,,(k), and Au, is the perturbation in the u, state. 
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The implementation of artificial neural networks trained via the back-propagation 
algorithm is described next. The implementation makes use of the methods described in 


this chapter, on feed-forward and recurrent networks structure, described in Chapter II. 
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IV. NEURAL NETWORK REPRESENTATION OF 
DYNAMICAL SYSTEMS 


It is interesting to examine how neural networks can interact with, act upon, or 
represent dynamical systems and processes. It seems evident that neural networks can 
have an important role in dealing with dynamical systems. There are already numerous 
examples of the application of neural networks to dynamical system prob'ems, and some 
examples can be seen in references [MSW90], [NP90], [NP91], [SS89], [HWB92], 
[SM92]. This research deals with the emulation and control of non-linear dynamical 
systems or processes using artificial neural networks. In this chapter a class of non- 
linear multi-input multi-output dynamical systems is developed. The affiliation of such 
systems with neural networks is then presented, followed by the special case of a linear 
time-invariant system. It is shown that a single network can represent a multitude of 
systems from the same class. A brief survey of previous work in related topics is 
presented, as well as motivations for extensions in this work. 

Narendra and Parthasarathy have presented extensive work, [NP9Q] and [NP91], 
on dynamical systems and their neural-network representation. They dealt with Auto 
Regressive Moving Average with eXtemal input (ARMAX) model structures and mainly 
with single-input single-output (SISO) systems. Based on the models described in 
references [NP90] and [MSW90, pp. 126-128], another presentation is suggested for 


multi-input multi-output (MIMO) state-space models. 
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It is possible to transform ARMAX form to a state-space form, since state-space 
representation is preferable. State-space representation is especially useful in giving 
insights into the physical mechanisms of the system, which can usually be more easily 
incorporated into state-space models than into ARMAX models [LL87, p. 81]. Another 
reason that working in state-space became very popular, e.g., in modem control theory, 
is that state-space models are more reliable in computer analysis [GLLT90, pp. 
1-67, 1-68], since the numerical algorithms are more stable. In fact, for MIMO systems, 
the state-space representation is the only convenient model. 

Most physical systems are more easily described in the continuous time domain, 
simply because most laws of physics are thus expressed. However, in many real-life 
implementations, a discrete representation is used and the modeling developed here will 
be based on a discrete representation. A similar description for the continuous time 1s 


achieved by replacing the difference equation by differential equations. 


A. CHARACTERIZATION OF NONLINEAR DYNAMICAL SYSTEMS 

Four classes of unknown nonlinear discrete dynamical systems were suggested by 
Narendra and Parthasarathy [NP90]. The models are extended here for MIMO 
dynamical systems with the addition of an output measurement equation. Common to 
all of them is that the state at t=(kK+1) is derived directly from the state and excitation 
functions at time t=k. If previous states are needed to describe the system, the prior 
states are augmented, forming a new state vector. The models are assumed to have the 


measured output of the system as a linear function of the state vector. 
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The nomenclature adopted here uses symbols that are usually used in control 
theory. In addition, upper-case italic bold letters followed by brackets represent vectorial 
functions, subscripted with / for a linear function, and with ai for a non-linear function. 
Upper-case bold letters with or without succeeding parentheses, represent matrices. The 
state vector is denoted as x, the measured output vector is denoted as_y, the input 
sequence is u, the discrete time is represented by k, and @ are the unknown parameters 


of the system. 
1. Model I 


x (k+l) = F,[k,O)x® + G,,[k0,,u] 5 x(0) = x, 
(4.1) 


y(k) = C(k) x(k) + D(k) u(k) 


In this nonlinear model, the state vector depends linearly on its prior state and 
nonlinearly on the input signals. The output is a linear combination of the state and input 
vectors. This model is very common in many mechanical systems where the inputs are 
subjected to nonlinearities such as saturation, backlash (dead zone), and switching 
(bang-bang). 

2. Model 


X(k+1) = F, [K,0,,x (k)] + G,[k,@ Ju(k) ; x) = x, 
(4.2) 


y(k) = C(k) x(k) + Dik) u(k) 


Here the nonlinear function F., acts on the previous states, while the inputs 


are linear. This is a very popular model in systems where kinematics are involved 
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(nonlinear trigonometric functions), systems with viscous friction, e.g., airplanes. This 
model or its linearized version, is very common in implementations of control actions, 


u, on the system F. 
3. Model Ii 


x (k+l) = F,,[k,6,,x(0] + Gy[k,0,,u@] ; x) = x, 
(4.3) 


y (kK) = CA) x(k) + Dk) ul’) 


In model III, both the input values and the previous state vector are subjected 
to two separate but additive nonlinear functions. Real-world examples for this model can 


be simply a combination of examples from models I and I. 
4. Model IV 


x(k+1) = H[k,0,,x(K),u(Q]; xO) =x, 
(4.4) 


y (k) = C(k) x(k) + Dik) u(k) 


This is the most general form that represents the current state as a nonlinear 
function of both the previous states and inputs combined. Although it is the most 
general, it is the most difficult to deal with, and tn practice, approaches described in 
models I to HI are mostly used. 

For every suggested model there is a corresponding neural network structure. The 
Structure characterizes the physical behavior of the model, and is trainable by the 
generalized delta rule. The representations and their generalizations are described in the 


following section. 
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B. NEURAL NETWORK REPRESENTATION OF NON-LINEAR 

DYNAMICAL SYSTEMS 

To represent linear or non-linear dynamical systems, recurrent networks must be 
used. The issue of building a proper recurrent network structure to represent, emulate 
or identify a non-linear dynamical systems has been treated by researchers in different 
ways. Narendra and Parthasarathy worked with SISO systems [NP91]. Werbos [WP90] 
and Gherrity [GM89] discuss training aspects for a recurrent networks. Seidi and Lorenz 
suggested a general structure by which a recurrent neural network could approximate a 
non-linear dynamic system, from the general form of model IV, where the measurement 
equation was a non-linear function of the states only [SL91]. Unfortunately, they did not 
address the problem of finding a weight adjustment algonthm to insure convergence of 
the system. 

A structure containing hidden recurrent layer, was also investigated by the author, 
in terms of training via back-propagation method. In the general case, when the output 
does not have a full presentation of the state vector, successful training via the back- 
propagation method is not guaranteed. The reason for the inability to learn to a 
reasonably small error is that during training the values of the state vector are unknown. 

The research conducted here classifies dynamical systems in two categories. First 
by the order of the input and output vectors, and the second by the relation of the output 
vector to the state vector. Each category 1s then subdivided as follows: 

Category A: Input-Output Orders 


1) Single-input, single-output (SISO). 
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2) Single-input, multi-output (SIMO). 
3) Multi-input, single-output (MISO). 
4) Miulti-input, multi-output (MIMO). 
Category B: Output to State vectors relation 
1) Output vector contains full state representation. 
2) Output vector contains partial state representation. 

For emulation, the dynamical system has to be completely reconstructible (also 
defined as constructible in the literature) in the linear case [KS72, pp.65-79,462-465], 
or realizable in the non-linear case [IS89, pp.129-144]. An examination shows that 
systems belonging to category B.1, i.e., systems containing full state output are rather 
trivial cases to emulate or identify using artificial neural networks. An example of linear 
time-invariant system with full state output is demonstrated and analyzed later in section 
IV.C. Systems from categories A.1 and A.2, 1.e., single input systems, are also rather 
straight-forward, since there is a theoretical solution for the one step ahead prediction and 
control [LL87, pp.51-65], that is easily adapted to neural network implementation. Most 
examples in the literature deal with models from or converted to categories A.1, A.2 and 
Bl. 

The main thrust of this work was to find a network structure to represent MIMO 
dynamical systems, trainable via back-propagation algonthm. The general form in which 
an ANN is trained to emulate or identify a plant is presented in Figure 4.1, where the 
inputs and desired values are generated, either by simulation or experimental data, and 


presented to the network in parallel. The network is trained until the error at the output 
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reaches an acceptable level. The specific inner structure of the network is now 


described. 





Figure 4.1: A Neural Network Emulation of Dynamical System 


In constructing an ANN it is important to take into consideration the structural and 
trainability aspects, which have to comply with theory. To better explain the structure 
and operation of an ANN used to emulate a dynamical system, first the solution is 
introduced followed by the incentives that motivated the approach. A representative 
network used to emulate dynamical systems represented by models I-Hi as described in 
section IV.A is given in Figure 4.2, where the circles denote processing elements, solid 
lines connecting PEs are adjustable weights, and the dash-double-dot line (—- -—) 
represents a fixed unity-valued weight. The structure is for any MIMO system. A 
system with three inputs and two outputs was selected for demonstration purposes only. 
The configuration retains the "black-box" form of the emulated system, by interfacing 


with the environment only with current values. 


1. Structural Aspects 

The structural considerations are considered in sequence. Systems from 
category A.4 - MIMO are considered, noting that categories A.1-A.3 are subsets of the 
former. The discussion starts with systems from category B.1 - full state output, such 
that for the measurement matrix C, rank (C)=n,, and n, 2 n,; where n, and n, are the 
dimensions of the output and state vectors respectively. For this case it 1s possible to 

talk in terms of y instead of x. For example, model II can be written as 
y (k+l) = F,,[k,8,,y(K)] + G,,[k,9,, 40] (4.5) 
where F’, G relates to F and G as y of eq. (4.5) relate to x in eq. (4.3), given that x(k) 
is the unique solution of C(k) x(k)=y(k)-D(k) u(k). The unknown parameters of the 


functions F, G are denoted 8, and 0, respectively. 

Each of the linear and nonlinear functions presented i models I-IV can be 
represented by a multi-layer neural network. The functions F,, F, are mapped by a 
network from the class Ne datacs n»,- Similarly, the functions G,, G,, are mapped by 
ae where 7, is the dimension of the input vector, and H, is mapped by 
ING . The network combined from feed-forward and the recurrent parts for u 


(n +n) Phy pene sft 


and y respectively. In this case, 
y(k+l) = N_[y(’)] + N,[u(’)] (4.6) 
where N is a recurrent network shown on the right of Figure 4.2, while N, is a forward 


network as depicted in the same figure on the left. The network N. is an extended form 


of the Hopfield recurrent network represented in a block diagram in Figure 2.5, where 


z' denotes a one-step delay. Note that Hopfield network is from the class N', ,» where 
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the input is fully connected to the output layer and the output values are fed back to the 
input layer via a one-step delay line. 

The more difficult case is for systems from category B.2, where the state 
variables cannot be computed from the output vector. For systems from category B.2, 
The measurement matrix C has rank (C)<n,, making the determination of N. more 
difficult. From linear systems theory it is known that the system needs to be 
constructible [KT80, pp.90-103], or reconstructible [KS72, p. 67] in order to be able to 
predict succeeding outputs. In non-linear systems, realizability is required [IA89, 
pp.131-144]. Reconstructibility means that present output values can be predicted from 
past values of the output. The need for previous information leads to introducing of 
regression layers. 

The inherent structure of the regression layer acts as a delay-line to provide the 
network with previous measurements of the input and output of the system. The order 
of a regression layer is defined as the number of previous measurements taken. The 
order of a regression layer is determined by the system which the network has to 
emulate, and sometimes may be assigned theoretically (see also minimal realization in 
Chapter VI.D). Examples of two regression layers of order three are shown in Figure 
4.2. On the left side the data from the input layer is regressed, and on the right side of 
Figure 4.2 the output values are regressed structuring a recurrent network. The 
regression layers become the practical input layer to the network, even though they are 
represented as an internal layers. The delayed values can be viewed as additional inputs 
to the network, to help create a proper relationship between the input and the output of 


the network. 
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The network is structured to correspond with the model chosen to represent 
the system. In the case of models I-III, represented by Figure 4.2, the effect of the 
input, described by G, is represented by JN; in the left branch of the figure, and the past 
output influence, described by F is represented by N at the mght branch. The 
resemblance to the physical behavior gives more insight to the system when analyzing 
the trained network and the weight matrices. Another advantage is gained when using 
the emulated system for training a peripheral network, such as training a controller 
network for the system. In the given structure, the back-propagated error that appears 
at the input of the emulated system is due to the input’s effectiveness, as opposed to 
errors spread over a large number of nodes, part of which are not relevant to the 
system’s inputs, as happens when one uses the generalized-delta-rule. From numerous 
simulations it was found that a network structure similar to the mathematical model 


yielded better results. 


, Regression —_ | 
Layers 





: oo N pa | 


Figure 4.2: Artificial Neural Network Emulator for Models I-U1 
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A suitable network configuration was established based on physical behavior 
of the system, and theoretical considerations for MIMO dynamical system. A suitable 
training procedure is next discussed. 

2.  Trainability Aspects 

It was shown in Chapter I.E that an ANN used to approximate a 
function/process, uses a squashing activation function in the PEs. Ifa squashing function 
is used, then the network has bounded output. In a physical system it is usually 
desirable, and practical also to have bounded inputs. Therefore, a dynamical system that 
is to be represented by an ANN should normally be a bounded-input bounded-output 
(BIBO) system over a time interval [7,2]. A system which is BIBO over time interval 
is denoted as bounded system ( see definition 3-36 in [AF66, p.101]). Note that a 
bounded system does not necessarily impose BIBO stability. In fact, unstable linear 
systems do belong to bounded systems, and can be emulated, as shown by Bertrand and 
Collins in [BC92] and [BD91]. It is worth mentioning that using a squashing activation 
function is not necessary. Indeed, to represent linear systems, better results are often 
obtained using a linear activation function [SR89], [BD91]. 

For bounded systems, two training schemes are possible, parallel and series- 


parallel schemes [NP90, pp.13-14]. For example, assume the dynamical system 
y*(k+1) = O[y%(k),u(k)], where Q is a vectorial bounded function to be emulated by N, 


with y(k+1) as the output of the network. The two schemes are herein described, using 


this example. 


48 


In the parallel-scheme the current network output is transferred to the output 
regression layer. The intermediate network approximations are used while training and 
in the recall phase, resulting in a process of the form 

y(k+l) = N[y®,u@® = 
The parallel-scheme has not been shown to be globally stable, even in the linear case 
[MSW90, p. 130]. 

The series-parallel scheme, uses the desired / known values, y’, in the output- 

regression layer while training, making a simple feedforward network. During the 


training the procedure takes the form of 
y(k+l) = N[y4(,u@) (4.9) 


Since the input and output values are assumed to be bounded, and there is no feedback 
loop containing non-linear elements, the standard back-propagation algorithm can be 
used. After training is complete y(k) ~y°(k), the parallel scheme is used in the recall 
phase, i1.e., the network outputs are regressed. The sSeries-parallel scheme has been 
shown to be stable [MSW90, p.130], and thus it will be used in this research. 

This section established a suitable structure and training scheme for an ANN 
emulation of a non-linear MIMO dynamical system. A simple and unique result arises 
in the linear time-invariant systems. The power of ANN in identification and emulation 


of linear systems is demonstrated in the following section. 
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C. ARTIFICIAL NEURAL NETWORK REPRESENTATION OF LINEAR 
TIME-INVARIANT SYSTEMS 
The linear time-invariant (LTI) system is a special case of either model I, II or I, 


where both F and G are linear functions. Consider the dynamical system 


x(k+1) = Fx(k) + Gu(k); x(0) =x, 
(4.10) 


y (kK) = Cx(k) 


where C is full rank and invertible when the number of outputs is equal the number of 


state variables. That specific system can be represented by a fully connected neural 
network from the class N’, , with linear activation function for the PEs, where the 


output layer is fully connected to itself, as seen in Figure 4.3 for an example of a system 
with four inputs and three outputs. 

For any network it is unknown in advance to which values the weights will 
converge to, and whether the solution is unique. The final results strongly depend on 
the structure of the network and the initial weight randomization. In the LTI case 
however, not only is convergence guaranteed but it is also unique; perhaps surprisingly, 
the weight matrices have a simple physical meaning. Examining the weights after the 
network converged, gives insight to the behavior of the system. In fact, if the system 
is known the weights are known a prion, and no training is needed. Using the fact that 


C is invertible, eq. (4.10) can be written in terms of y only as 


y(k+l) = CFC y(k) + CGu(k); y(0) = Cx, (4.11) 
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In several tests with LTI systems, the weights converged as expected to W''!=CG, and 
W"I=CFC', with residual error of zero if numerical rounding errors of a computer 
simulation were ignored. 
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Figure 4.3: A Neural Network Representation of LTI System 


D. EMULATION OF MULTIPLE MIMO SYSTEMS IN A SINGLE ANN 

Throughout the chapter it was shown that artificial neural networks have the 
capability of approximating non-linear dynamical systems. A new method is now 
proposed to extended the capability to emulate multiple systems from the same class in 
a single network. 

The different systems are assumed to have the properties discussed in subsection 
B. The dynamical systems belong to the same class and have the same input and output 
order. Additional type nodes, as many as the different models, are introduced at the 
input layer. At the presentation of information of a specific model, the corresponding 


node takes the value of one while all the rest of the type nodes are valued minus one. 


S1 


The training is in series-parallel scheme. To ensure proper training of all systems 
involved, the following conditions were obtained: 
1) Random presentation of each system. 
2) Balanced presentation of all systems. 
3) Persistent excitation in the domain of interest. 

The domain of interest in dynamical systems refers to both frequency and 
amplitude, over which the systems are to be mapped or emulated. The subject of 
persistent excitation will be discussed in more detail in Chapter VI.B, and is important 
for the single system as well. 

The implementation of the three conditions in this research for different dynamical 
models was as follows: 

@ Model is changed every presentation. 
@ Model number is selected for a random sequence without repetitions. 
@ Input signal is random-binary-sequence. 

The input signal is band limited to the desired spectrum, and amplitude of +1, with 
corresponding scaling of the models. 

Until now the discussion concentrated on emulation of dynamical systems by 
artificial neural networks. The next chapter discuss the role of ANN in adaptive control 


of dynamical systems. 
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V. ARTIFICIAL NEURAL NETWORKS IN ADAPTIVE CONTROL 


The control of LTI systems is a mature field with well established foundations on 
conditions and techniques, which define how LTI system can be controlled to yield 
desired responses. On the other hand, requirements for stable controllers of non-linear 
and time-varying systems are not generally available. In the last forty years, there has 
been significant progress in developing methods for identification and adaptive control 
of linear systems with unknown or time-varying parameters. Adaptive control techniques 
were based on the assumption of reasonably slow-varying changes in the parameters of 
the system, and depend on methods from linear estimation and control theory [AW89]. 
A recent trend of detection abrupt changes [AD91], [DS89], may be introduced to the 
adaptive methods for more rapid changing systems, but still staying at the linear regime. 
In the last decade, the nonlinear properties of ANN have been used to address 
identification and control of non-linear systems as indicated in the introduction. 

In this chapter a brief description of classical adaptive control schemes is given, 
and a description is also given of how the neural networks employ the classical schemes. 
A new method to form a neural-network adaptive controller is suggested. The idea 1s 


based on multiple-model approach together with the generalizing capability of ANN. 


A. BASIC CONCEPTS OF ADAPTIVE CONTROL 
The general idea in adaptive control is to find a regulator that will compensate for 


uncertainties, small or slowly changes in the controlled plant to form an overall stable 
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system. There are five basic approaches to adaptive control as described by Astrém and 
Wittenmark [AW89] including: 

1) Constant gain robust control. 

2) Self oscillating adaptive controller. 

3) Gain scheduling. 

4) Self tuning. 

5) Model-reference adaptive systems. 
Although the model-reference method is used in this thesis, a brief description is given 
of the other methods. 

Robust control is based on a linear state feedback for linear MIMO systems to form 
a large stability margin, good disturbance attedtialon and low sensitivity to small 
variations in the plant. This field is covered by the H, and the frequency weighted LQG 
theories. Robust controllers designed by means of H,, or H, schemes have been applied 
in many areas of multidimensional systems [CS88] as well as aircraft [HG90]. 

The self oscillating method is also a feedback based mechanism, utilizing a relay 
in the feedback loop. The bandwidth is automatically adjusted to be as high as possible 
and to a gain margin of approximately two. However limit-cycle oscillations occur. 
This bang-bang type of controller is an optimal controller for some systems, in other - 
human operated systems, it may be unacceptable. 

Gain scheduling is basically an open loop approach. The gains in a feedback loop 


are changed according to some measurable parameters in a pre-designed classical 
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controller. This method was originated for flight control systems to compensate for 
Mach number and dynamic pressure variations. 

The self-tuning regulator is based on a two-step procedure. Initially, there is an 
estimation process for the parameters of the plant. Based on the estimated parameters, 
the controller parameters are then updated under some design law. The parameter 
estimation is usually done recursively, using the results as if they were true. The 
controller design law as well as the estimation procedure vary ard are problem 


dependent, although should tightly correspond to the plant structure. 


1. Model-Reference Adaptive Control 

Model-reference adaptive controllers, sometimes referred to as model- 
following, involve methods to adjust the parameters of the controller so that the overall 
closed loop system will behave closely to a prescribed model. It is not always possible 
to achieve perfect following, but the overall stability properties are determined by the 
reference model. Caution must be taken to ensure a stable parameter adaptation 
procedure. The adaptation is based on some norm of the difference between the plant 
output and the reference model output. 

The model-reference, as well as self-tuning approaches can be further 
subdivided into direct and indirect adaptive controllers. In direct control the parameters 
of the controller are changed directly based on some measure of the output error, as 
shown in Figure 5.1. Alternatively in indirect control first the parameters of the 


unknown system are estimated, and are used to update the parameters of the controller 
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as presented in Figure 5.2. Both direct and indirect approaches result in a non-linear 


controller even for LTI systems. 


DESIGN OBJECTIVE 
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Figure 5.1: Direct Adaptive Control 
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Figure 5.2: Indirect Adaptive control 


To apply adaptive control techniques on LTI systems with unknown 
parameters, prior information should be available on the plant transfer function, such as 
the order and relative degree, and whether the zeros lie within the unit circle. The 


information is needed to determine the existence of the desired controller as well as the 


56 


reference model. It is reasonable to assume that for non-linear system even more 
information might be needed for theoretical treatment. Artificial neural networks can 
then be of help to find a suitable controller for the non-linear systems, as presented in 


the following sections. 


B. ADAPTIVE CONTROL USING ARTIFICIAL NEURAL NETWORKS 
Artificial neural networks are used to control dynamical systems in a variety of 
ways. Whether it is an imitation of classical methods, or utilizing the capability of ANN 
to represent a non-linear relationship between input and output variables even in the case 
where analytical depiction does not exist. An approach such as inverse dynamics control, 
which is not very practical in classical methods, is quite easily done using neural 
networks. In any case, new ways of application are opened to facilitate control 
procedures that otherwise would have been very tedious, if at all possible. This section 


summarizes different ways of using neural networks for control. 


1. Imitating Classical Methods 
Many of the applications of neural networks controllers are basically an 
imitation of classical methods. Robust controllers derived using H,, method were 
emulated using neural networks by Bertrand to control an unstable LTI X29 aircraft 
model {[BD91]. Neural networks were also used in the gain-scheduling approach as done 
by Burgin and Schnetzler for X15 pitch control [BS90], and for selecting proper gains 
and time-constants in a linear controller of the F/A-18 pitch controller by DiGirolamo 


and Donley [SM92]. Nguyen and Widrow used an ANN for the truck backer-upper 
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problem in a self-tuning fashion [NW90]. Ungar addressed the problem of bioreactor 
control using neural networks in term of self-tuning regulator [MSW90, pp.387-402], and 
Chen directly implemented ANN for self tuning adaptive controller for non-linear SISO 
plant [CF90]. The most common approach is the model-reference controller, since it has 
the least restrictions on the knowledge of the controlled plant. Examples of model- 
follower neural networks controllers for non-linear dynamical systems have been given 
among others by Narendra and Parthasarathy [NP90], and by Saerens and Soquet [SS89] 
for a variety of different plants. Comparisons between neural controllers and classical 
adaptive controllers [SSRB90], [KC90] show that while it is worthwhile using neural 
networks for control, there are still areas to investigate in utilizing neural controllers in 
practical and theoretical issues. The model-reference approach adopted in this research 
is discussed later in more detail. 

The method of inverse-dynamics control is a quite complex procedure as 
Classical approach and is rarely used [IA89, pp.181-183 and 234-243]. Neural network 


implementation of inverse control is more natural, as discussed next. 


2.  Inverse-dynamics Control 
Inverse-dynamics control is based on supervised learning, where the inputs 
to the ANN are the outputs of a system, and the desired outputs are the inputs to the 
system. A schematic of the procedure is given in Figure 5.3.a. When the network is 
trained and connected in series as in Figure 5.3.b, exciting with the desired output 
response, y*, will produce the corresponding control action such that the output of the 


system, y, is Close to the desired output value, i.e., y = y’. 
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The inverse plant was used by Bertrand to control a LTI unstable X29 aircraft 
[BD91], and by Teo for a non-linear SISO systems [TC91]. Kawato used a slightly 
different scheme of inverse dynamics for control of a multi-joint arm trajectory with a 


procedure to update the inverse plant when the dynamics of the system changes [MSW90, 


pp. 197-228], and [KIS88]. 






a) training 


INVERSE 


PLANT 





b) controlling 


Figure 5.3: Direct Inverse Control 


There are two major drawback of the inverse modelling. The first assumes 


that an inverse-plant exists, which imposes a limitation on the family of systems that can 


be controlled using this method. Second, it is rather sensitive to changes in the original 


plant, and on-line adaptation is rather long. Kawato also addressed the problem when 
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the inverse was not unique, 1.e., there is nO One-to-one mapping, as when different inputs 
produce the same output. The non unique mapping arises also when there are more 
actuators than degrees of freedom. But, when such an inverse mapping does exist, it 1s 
a quite straight forward and promising procedure. Indeed the H,, theory often involve 


elements of an inverse of the plant. 


3. | Model-Reference Control Using Neural Networks 
Model-reference adaptive control is based on reference model that specifies 
the behavior of the overall system. The controller parameters are updated based on the 
error between the outputs of the reference model, y’, and the plant, y. Figure 5.4 shows 
the general architecture of model-reference adaptive controller, where, r, is the command 
signal fed into the reference and the controller, and u, is the control actions obtained 
from the controller and excite the plant. Two loops are involved. The inner loop 


specifies the feedback controller, whereas the outer loop adapts the controller’s 


parameters. 
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Figure 5.4: Model-Reference Adaptive Control Architecture 


The model-reference approach overcomes some of the limitations of the 
inverse plant control. First it does not depend on the existence of an inverse plant. 
From the beginning of the training the control actions are in the desired trend to bring 
the plant to follow the reference model, 1.e., there is no need for a preliminary off-line 
learning phase. The procedure is directed to work on the desired domain which is the 
output of the plant, despite the fact that the controller does not have direct access to plant 
output, which is usually needed for training artificial neural networks. Since it is 
continuously learning to supply control actions in the desired direction, it can be used for 
varying systems. The problem is to ensure stable learning procedure over the possible 
changes in the system, which is not guaranteed. The adaptation to the new plant is also 
a long, time-consuming process, in particular when stable learning is embraced, thus 
questioning the capability of real-time application. The procedure suggested in this work 
addresses the stability and real-time limitations as discussed in section V.C. 

Model-following can be applied in neural networks both in the direct and 
indirect schemes. The indirect model reference control scheme as presented in Figure 
5.5 involves a preceding or interlacing phase of emulation of the plant. First the 
emulation network is trained and the weights are frozen. At the second phase the error 
at the output with respect to the reference model, is back propagated via the emulated 
plant and through the controller network, updating the weights of the controller to drive 
the output of the plant to follow the reference model. This approach was used by Psaltis, 
Sidris and Yamamura for coordinate transformation [PSY88], and by Jordan for a six- 


degree-of-freedom manipulator [JI89]. 
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Figure 5.5: Indirect Model-Reference Adaptive Control 


Direct model-reference adaptive control overcomes the need to first emulate 
the plant by utilizing the back-propagation through the plant described in Chapter I.E. 
A schematic of the direct model-reference controller is given in Figure 5.6. The neural 
network controller 1s trained based on an error, transformed backwards through the plant, 
and effects the plant output directly. This method was used by Psaltis, Sidris and 
Yamamura for coordinate transformation [PSY88], and by Saerens and Soquet for 
different constant dynamical systems [SS89]. Ha, Wei, and Bessolo supposedly used the 
method for control of a linearized longitudinal approximation of the F16 aircraft. 

The direct model-reference adaptive control with the extensions described in 
section III.E is the method used in this research. The specific implementation of the 


control scheme is presented in the coming section. 
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Figure 5.6: Direct Model-Reference Adaptive Control 


C. ANN ADAPTIVE CONTROLLER FOR TIME VARYING DYNAMICAL 

SYSTEMS 

The ultimate function desired from an adaptive controller is the ability to control 
a dynamical system under changing environmental conditions as well as changes in the 
plant itself. The dynamical systems are assumed to be output-controllable under all the 
changed conditions for the purpose of this design [OK90, p.706]. The adaptive control 
scheme is of the model-reference type. The reference model represents the desired 
response of the overall controlled dynamical system. The reference model is usually 
described by a linear dynamical system from r€R! to y“ER°, simply because the 
properties of linear systems are well defined. It is assumed that the reference model 1s 
chosen such that for a given bounded command signal, r, there exists a control action, 


u€v, exciting the plant, such that the output y follows y*. The control actions spanned 
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by UER™ and are bounded in accordance and as required by the plant. This in tum 
makes the controller be a bounded system by itself. Under the specified control actions, 
u€¥V, the neural network controller treats the "non-exact" tracking, by reaching the best 
controller in a least-mean-square sense. The changing plant is also assumed to remain 
within the same class of dynamical systems and retain the same inputs and outputs, 1.e., 


/,m and n are constant for the problem. 


1. Design Philosophy 
The ANN based controller design is adopted from classical control. A 
feedback controller N,,€ R°->R™ and a forward shaping filter N,€R'>R™, as shown in 
Figure 5.7 schematically. The two parts of the controller can be combined into one 
network which accomplishes both functions, having as inputs the command signal r, and 
the plant output y, resulting in the control actions u. The dash-dot (---) line in Figure 
5.7 indicates the combined controller, denoted N.. The training of the network is done 


as described in Figure 5.6 for the model reference approach. 
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The controller in the model-reference scheme is continuously adjusting to 


change the parameters in the network to drive the plant to follow the reference output. 


As mentioned earlier, it is desirable to shorten the transition time of the adaptation 
between changes of the plant. 

Training can be accomplished by defining a finite set of probable changes, 
each with its corresponding reference model and a resultant controller. Each model is 
trained off-line. As discussed in the multiple-systems emulator, a single network is used 
to incorporate all of the controllers creating a unified controller. Inputs to the controller 
network are used to relate the controller to the proper plant or new plant. For linear 
MIMO systems described in a state-space observer canonical form, controllability implies 
identifiability [LL87, pp.119-123]. A sensible choice of physical information to describe 
the plant, even in the non-linear case, is past measurements of input and output of the 
plant (see Figure 5.8). Having the information about the current plant in physical terms 
enables the network to generalize a suitable controller for the cases where the real plant 
differ somewhat from the plant the network was trained on. Consequently the network 
is valid for a variety of changes with a finite number of pre trained conditions, unless the 
given change is in some sense linearly independent of the models used to train the 


network. The structure of the unified controller is herein presented. 


2. ANN Structure for the Unified Controller 
The unified controller can be introduces in a network very similar in structure 
to the multi-system emulation network. Figure 5.8 shows the general block diagram of 
the training. One can see that the information on which a plant is currently controlled 
is passed to the reference as well, enabling a possibly different reference model for every 


plant. The network N., not only trained as a simple controller, but also makes a 
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synthesis of type information to generate proper control actions to the corresponding 


plant. 





REFERENCE MODEL 
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Figure 5.8: Unified Controller Training via BPTP in a Model-Reference Scheme 


In the case of dynamical systems, the information about the fype of the plant 
happens to be the same as needed for creating the controller itself. Previous 
measurements of the commands, output and controls are provided inside the network by 
using regression layers. Figure 5.9 depicts the classical controller as appears in Figure 


5.7 with additional type information combined into a single network as in Figure 5.8. 
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Figure 5.9: General Structure of ANN Unified Controller 
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The example given in Figure 5.9 shows a controller network with one 
command input, two control actions as an input to a plant with three outputs, all of which 
are fed into regression layers of order three. Each of the sub-networks, Nj, and N, in 
Figure 5.9, may have its own dynamics based on the inputs to that block. The blocks 
of the feedback controller and the forward shaping filter incorporates networks of the 
from appears in Figure 4.2. The type information, i.e., the regressed value of the inputs 
and outputs of the plant are supplied to both sub-networks. The z’ block represent a 
delay of one time unit, and is inherent in the structure of the network. The structure of 
the network is further generalized in section VI.D to suit a wider range of problems and 


arrangements. 
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VI. DESIGN CONSIDERATIONS 


Before starting to work with neural networks, whether for emulation or control, 


several design considerations need to be taken into consideration. There should be some 


general knowledge about the system such as its interface with the environment. The 


following issues need to be addressed: 


I) 
2) 
3) 
=) 
5) 
6) 
/) 


8) 


Defining measurable input and output variables. 

Frequency and amplitude ranges of interest. 

Sampling rate. 

Exciting input character. 

Data record length. 

The class to which the dynamic system belongs, if possible. 
Normalization and scaling. 


Capabilities of the peripheral equipment on hand, such as measurement devices, 
analyzing tools - hardware and software, etc. 


According to the available information, a suitable network is constructed and properly 


excited. To verify that indeed the network has been trained to meet the design goal in 


the domain of interest, the following matters must be looked into: 


a) 


Selecting a design criteria. 


b) Define and execute validation issues. 


The above parameters can be called design variables, and the manner by which 


they influence the design of a neural network is discussed in this chapter. 
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A. DYNAMICAL SYSTEM PROPERTIES 

In order to have a successfully trained network, the dynamical system has to have 
certain properties, depending on the task performed, 1.e., emulation or control. The 
properties were defined earlier and are repeated here to make the discussion about design 
issues more comprehensive. The systems under consideration in this research are non- 
linear MIMO dynamical systems. For both emulation and control the systems are 


assumed to be bounded, not necessarily BIBO stable (see IV.B). 


1. System Properties for Emulation Design 
In order to successfully emulate a system, or multiple systems in a neural 
network the system has to be reachable in the linear case, or realizable for nonlinear 


systems (see IV.B). The measured variables are the exciting inputs and the plant outputs. 


2. System Properties for Controller Design 
In order to successfully control a system, or multiple systems, in a neural 
network the system has to be output-controllable (see V.B). For controller design the 


measurable variables are the command signals, the control actions and the plant outputs. 


B. SPECTRUM AND AMPLITUDE ISSUES 

In analyzing dynamical systems the frequency domain of interest has to be defined. 
For emulation, the spectrum is usually defined by the dynamical behavior of the plant. 
In control applications, the bandwidth of the reference model, or the spectrum over 
which the controlled system has to function, determines the frequency range. For non- 


linear systems it is important to define the amplitudes of the measured variables both for 
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emulation and control problems. The trained network will perform well within the 
specified amplitude range, but there is no guarantee on the performance outside this 
range. In addition to the above conditions, the sampling rate, exciting signal 
characteristics complying with the concept of persistency of excitation, both in frequency 
and amplitude, and consequently the amount of data collected need to be considered, and 
are discussed in Chapter VII while examining the test case. 

It is of great importance in neural networks that the inputs to the network 
have similar energy content over time for proper training of the weights. Since a 
network is also a non-linear device, matters of amplitude should be considered for proper 
operation in the test phase. Therefore, the input and desired output values are scaled and 
biased from the physical values to fall within the acceptable range of the network, usually 
[0,1] or [-1,1]. If the energy content of certain inputs differ from others, say greater, 
it should be further compensated by the scaling or by selecting a smaller learning rate, 
to let the input equally effect the network. The relative significance will eventually be 
seen by the final values of the weights. Failing to correctly scale and normalize the 
network in its interface with the environment, may prohibit successful training within a 


reasonable long training period, if at all. 


C. NETWORK STRUCTURE SELECTION 
The specific structure of ANN per problem is generally the choice of the designer. 
This selection is not unique as will be shown below. The structural constraints and 


options are discussed in this section. 
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One constraint is defined by the interface with the environment which specifies the 
number of input and output nodes in the network. The inner structure is more problem 
dependent. As explained in Chapter IV, it is desirable to have the inner structure similar 
to the class of model to which the system is suspected to belong. The resemblance gives 
a better insight to the system when analyzing the resultant network, 1.e., there is more 
to a network than just a mesh. In some cases it can lead to almost perfect representation 
with full physical interpretation of the weights, as demonstrated with the LTI systems 
(IV.C). In reality noise is present in the measurements, as well as in the process, in 
which case the "perfect" is in the mean-square sense. 

The similarity between the emulation and control network structures, led to the 
formation of a general structure, suitable to generalize the procedure for emulation, 
control and related processes. 

A general ANN structure is introduced in this section. The structure embodies the 
general interface needed for emulation, identification, classification or control of 
dynamical systems, given in state-space or time series form. The general structure is 
presented in Figure 6.1. 

The input layer may consist of any or all of the following parameters: 


@ Forced excitation - used to drive the dynamical system with the suitable inputs, 
discussed in subsection B. 


@ Unforced excitation - other inputs to the network, which are not the direct 
excitation, but may result from the forced excitation. For example, an output of 


a reference model. 


@ Type - information used to indicate a certain system. Used especially in multiple- 
system networks. 
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Any portion or combination of the input layer may be passed to the variable-order input- 
regression layer as shown in Figure 6.1. 
The output layer may consist of any or all of the following parameters: 


@ Dynamic output - used for varying outputs, such as the output of a dynamical 
system. 


@ Static output - used e.g., for classification, or other non-varying output parameters. 


Any portion or combination of the output layer may be passed to the variable-order 


output-regression layer. 
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Figure 6.1: A General Structure for Multi-Layer ANN 





The effector layer seen in Figure 6.1 is introduced if some intermediate values 
were needed to describe a more complex system. For example, combining two networks 
where one feeds the second, as in control problems. This layer is also regressible 


through the output regression or via the unforced excitation. 
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Different structural combinations are possible to solve the same problem. The 
inner structure, as discussed earlier is the choice of the designer, and problem dependant. 

A program to interface with NeuralWare’s software was written in C, based on this 
general framework for a multi-layer ANN, with easy accessible user-defined functions 
to accommodate a variety of systems, structural combinations and excitations. This 
USERIO program appears in Appendix D, along with a dedicated control strategy that 
monitors the execution of the USERIO program. 

This section established this structure of artificial neural networks in conjunction 
with non-linear dynamical systems. A new method is now proposed to define minimal 


realization of the inner structure. 


D. MINIMAL REALIZATION OF A MULTI-LAYER ANN 

Multi-layer networks consist of an input layer, an output layer, and at least one 
hidden layer. A typical feed-forward network with two hidden layers is presented in 
Figure 2.2. The layers are fully connected, and the connections form a weight matrix 
W. Theoretically, there is no limit on the number of hidden layers or the order of layer 
connection, 1.e., the connected layer does not have to be consecutive. A bias term is 
usually added to each PE with an associated non-zero weight that can be considered as 
a threshold of that PE. The number of inputs and outputs is dictated by the nature of the 
problem which the network has to solve or represent. The number of hidden layers, 
their arrangements, the number of PEs in each layer, and their connections are the 


designer’s choice. 
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Separate theorems express the existence of artificial neural networks as mapping 
devices from a finite dimensional space to another to any desired degree of accuracy with 
only one hidden layer, as discussed in section I.E. Kolmogorov’s theorem [HN90, pp. 
122-124] gives an upper limit to the number of hidden nodes needed for the mapping 
(2n+1, where n is the number of input nodes). However, for every mapping task a 
different unknown activation function is needed. It has also been shown [HSW89], that 
ANN are general approximators of a bounded function using an arbitrary squashing 
activation function, i.e., sigmoid or hyperbolic tangent, in an unknown but large enough 
number of PEs to get the desired accuracy. However, neither theorem really helps to 
build a “correct” structure of the network. 

An understanding of the physical nature of the problem may assist in choosing a 
“good” number of hidden layers and hidden nodes. Until the present, choosing network 
parameters was a kind of art usually done by tnal and error. A more rigorous way of 
selecting the number of hidden layers and how many nodes each layer should include is 
presented in this dissertation. 

Although the method cannot predict in advance a good choice for the number of 
layers and PEs, the rationale for going through the procedure of finding minimal 
realization of the network 1s to achieve the same accuracy with less network complexity. 
It is especially useful in recall mode by having less computation and less parameters in 
"the network, which leads to easier implementation and faster reaction. In training mode 
however, there is no distinctive conclusion on total time required for complete training. 


For off-line training there is no major impact but it may be useful in on-line processes 
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if it can be found that a much simpler network, which trains faster, can still do the job. 
Tests showed that usually the minimized network needed more iterations to converge but 
the overall training time has to be examined per task, and might even be shorter since 


less computation per iteration are involved. 


1. Selection of Number of PEs in Hidden Layers 

The method developed finds a minimal realization of the network, by means 
of number of hidden layers and number of PEs in each layer. It is assumed that there 
exist a network, trained to the desired accuracy, containing a large enough number of 
layers and PEs. The assumption is based on the results of Hornik, Stinchcombe and 
White [HSW89]. Thus, the network is assumed to have an adequate number of hidden 
layers, each containing large enough number of PEs, where every PE incorporates a 
squashing activation function. The idea is to find how many PEs carry most of the 
activation level from one layer to the next. It is done by examining the singular values 
of the weight matrix connecting the two layers. 

The singular values of a rank g matrix W€R™*", denoted by a,, are the 
non-negative square root of the eigenvalues of W'W ordered such that o, =0,2>... 209,. 


If q<n, (n-r) singular values are zero, 1.€., 04,;=0j4.=...0,=0. There exist two 


q 


unitary matrices UER™*™, VER" and a diagonal matrix L€ R™*" such that, 


w-urvt-|-7- (6.1) 
0 


where L,=dlag(o,,0,,...0,). This partition is called the singular value decomposition 
(SVD) of the matrix W. The greatest singular value o,,,, 1S given by 0. 

Since the weight matrix contains the information on how the data is passed 
from source layer to the target, the singular values of the matrix are indicative of the 
characteristics of the matrix and hence serve as a feature detection process. The specific 
features embedded in the weight matrix need to be extracted. To identify the features, 
a trace of the singular values of the weight matrix is calculated while training progress. 

Figure 6.2 shows a typical trace of the singular values of the weight matrix 
as a function of training iterations, where the unreduced target layer has seven PEs. The 
traces of each singular value is denoted by S1,S2,...87. The singular values carrying the 
important information will develop and can be identified from the rest that do not carry 
significant information. In a multi-layer network , the weights are usually given a 
random initialization within a specified range, carrying values other than zero. To be 
able to identify even very small but significant feature, the initial randomization should 
be small enough to let the small singular values emerge above the initial values. The 
trace marked by S6 is an example of a significant singular value that pops-up over 
singular values that initially were bigger but do not carry significant information. 

It should be noted that the singular values are traced continuously, and not 
necessarily sorted in descending order as defined by eq. (6.1). ANN that are trained 
using the back-propagation algorithm must have initial values other than zero. Small 
initial values are theoretically acceptable but may cause difficulties at starting the training 


process. 
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Based on the chante values an approach to select the number of PEs needed, 
is to define a threshold X, above which the contribution to the inputs to the target layer 
is considered significant, and set the number of PEs in the target layer as the number of 
singular values above the threshold. The threshold is usually related to go,,,, by a 
variation ratio - v as, 


Nop, = max {i|o,=A} (6.2) 


i=1,2,...,” 


PEs 


X= Lmax (6.3) 


v 


In the example given in Figure 6.2, the estimated number of PEs needed is four out of 
seven initially set. 

However, using the singular value decomposition to determine the number 
of PEs needed in the target layer, implies that initially the number of PEs in the source 
layer must be greater or equal to the number of PEs at the target layer. 

The plot of singular values can also serve as a useful tool to determine if the 
training is completed, by noticing that the singular values have reached their steady state 


values, and more training will not give any improvement. 


Singula: values 
fe) 





Figure 6.2: Singular Values of a Weight Matrix as a Function of Training Iteration 
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2. Selection of Number of Hidden Layers 
The number of hidden layers needed is determined experimentally in the 
following way. It is assumed that the network was previously constructed and 
successfully trained with an adequate number of layers. A layer is said to be redundant 
when the SVD plots of consecutive layers show the same number of significant singular 


values, and their traces are identical up to a scale factor. 


3. Selection of Regression Order 

The order of the regression layer, 1.e., the number of previous measurements 

taken into account as an input to the network, should first be based on physical or 
theoretical reasoning. If no information on the order were given, the initial regression 
layer would have to be of large enough order assumed to be sufficient, e.g., estimated 
order of the system. The number of input and output histories is determined by 
observing the weight matrix connecting the regression layer to the hidden layer. Figure 
6.3 shows a trained weight matrix from a regression layer of order eight, with three 
regressed values, to a hidden layer with five PEs. 
The height of each box is the value of the specific 


weight, and the same shading pattern refer to the 





same regressed variable. From Figure 6.3 it is 


1 


seen that the contribution of regressed values prior a 


+ R W 4 4 i W | 


to t-4 are negligible, and thus can be eliminated Figure 6.3: A Weight Matrix From 


leaving the regression layer of order four only. BS has aE Cen ahs 
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E. VALIDATION ISSUES 

To conclude the design considerations, one must define a design criterion, and 
supply the means to verify that indeed the goal has been achieved. It is not sufficient to 
require that some norm of the output error, for either emulation or control, will meet 
some threshold. It has to be proven that the task is fulfilled over the whole domain of 
interest. The concept of consistency of excitation holds here as well as some other 
measures to verify the results. 

In dynamical systems, a common way 1s to check the response in the time domain, 
to some known test cases and analyze them. Common time-domain tests are pulse, step, 
ramp or dipole commands, as well as checking the performance in presence of noise and 
evaluating the results statistically. A complimentary evaluation is done to verify the 
spectral properties. A frequency response is calculated in parametric or non-parametric 
means, and compared to previous knowledge on the system. The location of poles and 
zeros can be extracted, to give more physical interpretation on the system. 

With neural networks, another means to evaluate the results is done by examining 
the weight matrices of the trained network, assuming that there is some correspondence 
between the network structure and the physical behavior of the system. An example is 
given in the minimal realization procedure (VI.C), and by the LTI system emulation (see 
[V.C). For example, from the weight matrices one can estimate the number of previous 
measurements needed to represent the system, which gives information on the order of 


corresponding ARMAX model. 
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Most of the above criteria and techniques are used in this research to verify the 
results, as demonstrated in Chapter VII. 

Up until now, the theory behind neural networks and training mechanisms was 
established. It was shown how ANN can be utilized to emulate and control non-linear 
MIMO dynamical systems. Ways to construct a suitable network, train and evaluate the 
results were discussed. It is now appropriate to demonstrate all the above in an example, 
and show the strength and capabilities of artificial neural networks in emulation and 
adaptive control. The rest of the dissertation involves the presentation and the evaluation 


of a representative test case. 
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Vil. ANN EMULATION AND CONTROL DEMONSTRATION 


This chapter presents a test case to examine implementation of emulation and 
control using artificial neural networks as developed in the previous chapters. 

The F/A-18A fighter aircraft ae chosen as a demonstration MIMO dynamical 
system due to the availability of multiple control surfaces to a rather complex dynamics, 
see Figure 7.1 below. A damaged aircraft was assumed to represent a change in the 
system, which the neural network emulates and controls. Thus, the neural network 
forms a reconfigurable adaptive controller for an impaired aircraft. The closed-loop 
performance is dictated by a reference model that meets the conditions required by 
military specification MIL-F-8785C. The discussion continues with a description of the 
open-loop dynamic model of the F/A-18 followed by specific damage mechanism which 


were treated in this research. 


A. F/A-18A LONGITUDINAL DYNAMICAL MODEL 

The F/A-18A linearized longitudinal dynamics with asymmetrically-saturated 
control inputs, provides a 25" order dynamical system with six control inputs and three 
output variables. The dynamical model includes the airframe dynamics as well as 
actuator and sensor dynamics. The continuous state-space model was developed from 
data in [MDC82] for a flight condition of Mach 0.6 and altitude 10000 feet, in escort 


configuration. The aircraft equations of motion are taken from [MAG90, pp. 256-260, 


81 


296-298] in stability axes and trim conditions, assuming no gusts. The actuator and 
sensor models are based on the model developed by Rojek [RF386]. 

Figure 7.1 shows the F/A-18 with the control surfaces and describes the sign 
conventions for both deflections and aerodynamic coefficients. The terms in body axes 
are denoted (*), and (®), denote stability axes. The control surfaces affecting the 
longitudinal dynamics are the stabilators (dst), leading edge flaps (dle) and trailing edge 


flaps (dte). 
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Figure 7.1: F/A-18A Axes, Control Surfaces and Sign Conventions 
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The continuous state-space representation of the aircraft, which includes the 
sixteenth order actuator dynamics, fourth order airframe dynamics, and fifth order sensor 


dynamics is given in the 25" order augmented system as follows. 


x = Fx + Gu 
(7.1) 
y = Hx + Du 
The control vector, u, is the actuator deflection commands [deg], 
| {5,, b 6, b,. b.. 6, (7-2) 


where st, le and te denote the stabilators, leading edge and trailing edge, and the 
subscripts r,/ refer to the nght and left sides respectively. The deflections of the surfaces 
are limited to: 

@ Stabilators [+10.5°,-24°] 

@ Leading edge flaps [+34°  ,-3° ] 

@ Trailing edge flaps [+45° ,-8° ] 
The measurable output, y, consists of, 

y = {q n, a}! (7.3) 
where q is the pitch rate [deg/sec], the normal acceleration n, [g’s] and @ is the angle of 
attack [deg]. 

The numerical values of the aerodynamic stability derivatives together with the 
actuator and sensor dynamics appear in a MATLAB™ script file, f7SlongN.m, that 
calculates the augmented system, and are given in Appendix A, with the eigenvalues and 
the corresponding damping and frequencies appearing in Table 7.1. The program also 


calculates the matrices for the impaired model described next. 
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Table I EIGENVALUES, FREQUENCIES AND DAMPING RATIOS OF 
UNDAMAGED F/A-18A LONGITUDINAL CONTINUOUS MODEL 


S-Plane roots Damping Frequency Mode 


[rad/sec] 
-0.0054 + 0.05691 0.095 0.057 Phugoid 
-0.0054 — 0.05691 0.095 0.057 Phugoid 
-0.9786 + 2.63411 0.348 2.810 Short Period 
-0.9786 — 2.63411 0.348 2.810 Short Period 
-14.000 1.000 14.00 AoA Sensor 
-14.924 + 33.200: 0.410 36.40 Stab. Act. night 
-14.924 — 33.200i 0.410 36.40 Stab. Act. right 
-14.924 + 33.200i 0.410 36.40 Stab. Act. left 
-14.924 — 33.200i 0.410 36.40 Stab. Act. left 
-24.850 + 24.647i 0.710 35.00 TE Flap Act. nght 
-24.850 — 24.647i 0.710 35.00 TE Flap Act. right 
-24.850 + 24.647i 0.710 35.00 TE Flap Act. left 
-24.850 — 24.6471 0.710 35.00 TE Flap Act. left 
-26.900 1.000 26.90 LE Flap Act. night 
-26.900 1.000 26.90 LE Flap Act. left 
-62.127 + 85.020i 0.590 105.3 Stab. Act. night 
-62.127 — 85.020i 0.590 105.3 Stab. Act. night 
-62.127 + 85.020i 0.590 105.3 Stab. Act. left 
-62.127 — 85.020i 0.590 105.3 Stab. Act. left 
-82 .900 1.000 82.90 LE Flap Act. right 
-82 .900 1.000 82.90 LE Flap Acct. left 
-167.80 1.000 167.8 Pitch Rate Gyro Sensor 
-379.49 + 110.68 0.960 Ba) Normal Accel. Sensor 
-379.49 — 110.68: 0.960 Seiad, Normal Accel. Sensor 
-461.70 1.000 461.7 Pitch Rate Gyro Sensor 


B. DAMAGED F/A-18A MODEL 

The damage mechanism considered in this example is a malfunction in the actuator 
of the left stabilator. The left stabilator does not respond to actuation commands and 
remains fixed at trim condition. This is modeled by nullifying the corresponding actuator 


dynamics in the actuators model, as shown in Appendix A. 
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C. CONTROL REFERENCE MODEL 

The reference model for the model-following controller is based on the military 
specification for flying qualities of piloted airplanes, MIL-F-8785C. The general 
configurations is shown in Figure 7.2. The same reference model is selected for both 
damaged and undamaged F/A-18 models. 


REFERENCE MODEL 
(MIL-F-8785C) 
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Figure 7.2: Model-Reference Control Architecture for the F/A-18A 


The reference is a forth order linear model, that defines desirable behavior for the 
short period and phugoid modes [MS80, pp.12-16]. The frequency and damping are 
@ Short period - w,,=4.3 [rad/sec]; ¢,=0.6 
® Phugoid - w,,=0.02 [rad/sec]; ¢,=0.4 
The reference model is given in discrete state-space form as 


2(k+l) = F2(k) + Grtk) 
(7.4) 


yi(k) = H2(k) + Drik) 
The input command, r, is the pilot pitch stick commands p,. The output vector, y’, 


consists of, 


yt = {4 nf otf an 
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where q7 is the reference pitch rate [deg/sec], the desired normal acceleration n.’ [g’s] 
and a’ is the reference angle of attack [deg]. The reference model is scaled such that at 


full range pilot command p,=1 generates a normal acceleration of n,“=5g in magnitude. 


D. DESIGN PARAMETERS AND CONDITIONS 

This section gives the test case specific parameters and conditions for topics 
discussed in Chapter VI as design considerations. The fact that the basic aircraft model 
is linear makes it easier to apply the conditions set in Chapter VI, but neural networks 


have been shown to successfully work with non-linear functions as well. 


1. Plant Properties 
Both the undamaged and damaged aircraft models are bounded MIMO 
systems. The models are reachable as required for emulation, and output-controllable 


as required for control application. 


2. Spectrum and Amplitude Ranges 

The frequency band of interest is taken to be the same for emulation and 
control. It is chosen to cover the open-loop and closed-loop aerodynamic modes. The 
poles of the actuators and sensors are much faster and left out of this range. The 
frequency band of interest covers three decades from /,,,,=0.005 Hz to /,,,,=5 Hz. The 
sampling rate in this work is 50 Hz, which lets the Nyquist frequency be five times 
larger than the highest frequency of interest, and about ten times the short period mode 
and therefore aliasing is greatly suppressed. This sampling rate causes the discretized 


plant to become non-minimal phase. Non-minimum zeros are accounted for in the 
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control phase by resetting the simulation, to prevent modelling errors from growing 
unbounded. 

The exciting signal follows the concept of persistence of excitation. A 
random binary sequence (RBS) is used, shaped such that the energy content of the signal 
will be distributed across the desired spectrum. The exciting signal coincide with the 
sampling rate. Since the network is non-linear the exciting signal should cover the whole 
amplitude range of interest. The RBS signal was found suitable for persistently exciting 
the network range, and was set to span [-1,1]. 

The characteristics of the RBS signal are shown in Figure 7.3 in the time 
domain. Figure 7.4 shows the spectral characteristics of the RBS signal used for training 
(order of 500k presentations). The drop in the spectral energy above five Hertz was 
intended to limit the excitation to the desired range and prevent high frequency 
corruption. For testing the trained network it was impractical to use such long 
sequences, and the shortest sequence (16384) capable of spanning the desired range of 
[0.005 ,5] Hz was used. The spectrum of a testing sequence appears in Figure 7.5, which 
is less uniform then the training spectrum. It should be noted that it is very difficult to 
properly excite the very low end of the frequency range (f=0.005 Hz is a cycle time of 
200 sec. which is 10000 samples), therefore some mismatch might be expected in the low 
range due to insufficient excitation from practical reasons. 

The amplitude bound brings up the issue of scaling. The desired range of 
signals for a successful learning of the network is [-1,1]. Therefore the models of the 


plants and the reference model were scaled and biased from the physical range to the 
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network range. The biasing eliminated the non-linearity in the actuator’s deflection 
range. The non-linearity of the network is expressed in multiple models handled in a 


single network. 


Random Binary Sequence Sample. Ts#0.02 sec. 


Amplitude 
Sb 0 
ta 


oO 2 4 6 8 10 12 14 
Time [sec] 


Figure 7.3: Time Domain Sample of Random Binary Sequence 
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Figure 7.4: Power Spectrum of Training RBS Signal (512k presentations) 
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Figure 7.5: Power Spectrum of Test RBS Signal (16k presentations) 
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3. Network Structure for Emulation and Control of the F/A-18 Model 
The multiple-plant model belongs to the class of models which are 
represented by model III as of section IV.A. The corresponding structure was selected 
as appears in Figures 4.2 and 5.9 with the appropriate number of inputs and outputs, as 


will be seen in the result chapter next. 
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VUI. SIMULATIONS, RESULTS AND ANALYSIS 


This chapter describes the different experiments done with the aircraft models 
described in Chapter VII. The results are shown and analyzed, both in the frequency and 
time domains. The spectral plots of the neural network response and the desired output, 
used FFT analysis of the input and output sequences. The data was windowed using 
Hanning window. The system Bode plot is the analytic spectrum of the simulated 
system, which is used for generating the desired output. The neural network results are 
given for every relation between the input and output of the MIMO dynamical system, 
and compared with the theoretical and desired values. 

In the following plots, Plant #1 refers to the undamaged aircraft model, whereas 
Plant #2 refers to the damaged left stabilator actuator impaired aircraft model. The 
experiments start with the emulation process and are followed by the adaptive controller. 

The emulation results includes three parts. First the undamaged and damaged 
aircraft models are emulated separately in the proposed network structure that contains 
feedforward and recurrent sub-networks, corresponding to model III (see section IV.A 
and Figure 4.2). In the second part, a single network is used to emulate both aircraft 
models using a random balanced presentation procedure as described in section IV.D. 
The multiple-plant emulation in a single network is repeated with a minimized network. 
The minimal realization is done via the singular value analysis of the weight matrices and 
evaluating the minimal regression order by analyzing the weights from the regression 


layer to the first hidden layer as described in section VI.D. 
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Since the F/A-18 aircraft used in the analysis is a MIMO system with six inputs 
and three outputs, many plots are needed to completely describe the behavior of the 
system. Therefore the results shown in the following A, B and C sections are for 
representative variables, and the rest of the results, when appropriate, are gathered in 
Appendix B. The representative inputs are the right and left stabilators, noted in the 
graphs as dst, and dst,, respectively. These inputs were selected to emphasize the 
differences between the damaged and undamaged plants. Recall that the damage is a 
malfunction in the left stabilator actuator. The output variable selected is the normal 
acceleration n,. Note: the results due to the selected variables are not necessarily the best 


among the rest of the inputs and outputs. 


A. EMULATION OF A SINGLE MIMO DYNAMICAL SYSTEM 

The experiment started with an initial network that had regression layers of 25" 
order (same order as the plant, as might theoretically be required). Examining the 
weight matrix from the regression layer to the first hidden layer showed that the three 
most recent regressions carried most of the activation and the weights of prior 
measurements were rather small compared to the first three. Based on this observation 
the regression layers were reduced to be of order three, as in Figure 8.1 for the multi- 
plant emulation structure, and is the network used for the results presented here. The 
order reduction of the regression layer greatly reduced the complexity of the network and 
the training time, with no loss in the accuracy of the training. This was the first 


reduction of the network structure, but the hidden layers were not minimized. 
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The input regression layer denoted as ComReg, : a third order regression of the 
six control inputs containing 18 nodes. The input regression layer ts fully connected to 
the first hidden layer comhJ which is further fully connected to the second input hidden 
layer comh2, each containing 36 PEs. The second hidden layer is fully connected to an 
output-buffer layer, OurBuf, through a connecting Layer with unit-valued fixed 
corresponding weights, for visual separation of forward portion of the network from the 
recurrent portions of the network. In Figure 8.1, the connections tetween layers are 
symbolically illustrated by the connection from a single PE at a source layer fully 
connected to the target layer. Showing the full connectivity will make the drawing too 


complex to visualize. 
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Figure 8.1: Network Structure for a Multiple-Plant F/A-18 Emulation 
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The recurrent eroenctork Starts with a third order output regression layer OutReg 
which obtains the inputs via a feedback from the output layer. The use of the output- 
buffer layer was required due to limitation of the current software from NeuralWare that 
was used in this research. The OutReg layer consists of nine PEs and is fully connected 
to the first hidden layer of the recurrent portion outh1, which is further connected to the 
second hidden layer outh2, each containing 36 PEs. The second hidden layer is 
connected to the output through Layer, for better visualization. In the recurrent sub- 
network connections are also shown only schematically, since all the layers are fully 
connected. 

In this section each plant, the undamaged and damaged aircraft are emulated 


separately on two different networks structured as in Figure 8.1. 


1. Single Plant Emulation - Frequency Domain Results and discussion 

The examination starts with spectral analysis of the network compared to the 
desired model. Figure 8.2 shows the spectral response of the pitch rate q, the normal 
acceleration n, and the angle of attack AoA of the undamaged aircraft due to the nght 
stabilator dst, (upper plots) and the left stabilator dst, (lower plots). Each graph shows 
the analytic frequency response in the solid line, which is the Bode plot of the outputs 
due to the corresponding input, together with the desired outputs in the dashed line, and 
network output values in the dash-dot line, analyzed using FFT, while exciting the input 
with random binary sequence (RBS). Analogous plots for the other inputs to the system, 
i.e., the leading edge and trailing edge control surfaces, appear in Appendix B, Figures 


B.1 and B.2 correspondingly. 
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Thewticunee show almost perfect match between the system Bode plot, the 
desired output and network output frequency response over the range of 0.03 Hz to 5 Hz. 
The deviation between the network output and the desired output is mainly a problem of 
persistent excitation. Two factors contribute to that difference. First it is very difficult 
to properly excite at such low frequencies (0.01 Hz is a cycle time of 100 seconds and 
in 5000 learning iterations only one change is expected). The other factor 1s the inherent 
fact that as the frequency goes down, the number of presentations of data information due 
to that frequency goes down. The less a feature is presented to a network, the less 
trained the network is to correctly respond to similar excitation. 

The deviation of the network output and the desired output from the analytical 
Bode plot is due to the coarse resolution at the low frequency. Since the spectral 
responses of the network and the desired output are done using FFT analysis, the 
spectrum is linearly spread. In order to get better resolution over a range of three 
decades, the amount of data needed to be collected becomes prohibitive with respect to 
the capabilities of the current analyzing tools (software and computer memory). Due to 
the reasons mentioned above and the random nature of the input signal, at the low 
frequencies the general trend of the spectral behavior curve should be compared to the 
analytical system Bode, and not the exact numerical value. 

The two dominant modes of the aircraft at the frequency range of interest, 
i.e., the phugoid at the low frequencies and the short period at the high frequencies are 
clearly characterized by the neural network as can be seen in all the spectral plots in 


Figures 8.2 and B.1, B.2. 


94 


Figure 8.3 and the corresponding Figures B.3 and B.4 in Appendix B. 


presents the same information as Figures 8.2, B.! and B.2, but for the damaged aircraft 


trained in a separate neural network. The spectral results obtained for the undamaged 


aircraft hold for the damaged aircraft as well. In Figure 8.3 the damage effect is clearly 


seen in that there 1s no response of the aircraft to excitation due to the left stabilator 


(dstl) over the whole frequency range. The residual value of the network output at about 


-80dB in the lower part of Figure 8.3 is due to the initial weight randomization and 


numerical noise. The system Bode and the Desired value spectrum is minus infinity and 


thus do not appear on the graph. 


Looking at all six figures, Figures 8.2, 8.3 and Figures B.1-B.4 some general 


observations from the results are emerging. 


Ve 


The network learns the transfer function of each of the six inputs. This can be 
seen by comparing the upper and lower graphs of each figure. Since each figure 
shows the corresponding right and left surfaces, the results should be symmetrical 
for both, apart from Figure 8.3 where the damage is seen. The slight differences 
between the right and left surfaces of each of the three pairs indicate the 
individual treatment of each input. 


The network is able to resolve the similar input pairs for example left and nght 
Stabilators, leading or trailing edge flaps, although each is treated separately. 
This is seen by the very close resemblance of the spectrums of each pair, right 
and left inputs, which in the undamaged plant are identical, and also for the 
damaged plant except for the stabilators. 


The network is clearly able to identify the damage in the second plant as seen in 
the lower part of Figure 8.3. 


There is a very good agreement between the desired response and the network 


response over almost all the frequency range. The desired and network responses 
were calculated using the same random input sequence. 
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There is a very good agreement between the network response and the analytical 
Bode plot over about the upper two and a half decades. It can be said to have 
exact match over frequencies from 0.03 to 5 Hz. The differences are due to the 
random nature of the input signal used to generate the network response. The 
difference is also due to the difficulty to practically excite the very low 
frequencies, and the lower number of presentations of the phugoid mode relative 
to the short period mode. 


The networks have learned the spectral behavior of the two plants as presented 
in Figures 8.2, 8.3, B.1-B.4. 
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Figure 8.2: Spectral Response of Undamaged Aircraft Due to Stabilators 
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Figure 8.3: Spectral Response of Damaged Aircraft Due to Stabilators 
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2. Single Plant Emulation - Time Domain Results and Discussion 

The time domain analysis shows the response of the outputs of the system, 
i.e., g, n, and a, to a dipole test signal driving each input. Due to the multiplicity of 
inputs and outputs, and in order not to disrupt the flow of reading, representative figures 
follow here, and the rest are gathered in Appendix B. and referenced accordingly. The 
organization of all figures in the appendix is identical to the figures shown and described 
here. 

Figure 8.4 shows the normal acceleration response, n,, of the undamaged 
aircraft due to a dipole input command at the night stabilator dst, (upper plots) and at the 
left stabilator dst, (lower plots). The dipole presents a full deflection in the positive 
direction for one second followed by a full deflection command in the negative direction 
for a duration of another second, and then reset to zero deflection. The dotted line 
represents the dipole and 1s indicative of the switching times (the magnitude is +1 and 
out of scale). The solid line is the desired behavior, and the dashed line in the neural 
net output. The response of the normal acceleration due to leading edge flaps dle, and 
dle,, and due to trailing edge flaps dte, and dte, appear in Appendix B. Figures B.8-B.9 
respectively. The pitch rate response of the undamaged aircraft to a dipole command in 
all inputs is in Figures B.5-B.8. The corresponding angle-of-attack response appears in 
Figures B.10-B.12. 

Fach graph shows the desired and network output for the command. The 
error between the two is presented in the lower plot of each pair. Note that the output 


are on the order of tenths while the errors are at the order of 10°. The relative error of 
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all input-output relations is about 1%, which was the criterion for successful eacnine: 
It is interesting to see that the network which was trained by exciting all inputs 
simultaneously was equally trained on all inputs regardless of the individual behavior and 
amplitudes. Moreover, each input is interpreted individually as can be seen from the 
different error of similar control surfaces pairs (left and right), in Figures B.7, B.8 and 
B.9 (it is impossible to resolve this from the actual response because of the very small 
difference). 

The damaged aircraft was trained on a separate network. Figure 8.5 presents 
the normal acceleration response of the damaged aircraft due to the stabilators. The 
graphs of n, due to leading and trailing edge flaps appear in Figures B.16, B.17 
respectively. the corresponding plots of the damaged aircraft pitch rate response to a 
dipole command in the various control surfaces are in Appendix B, Figures B.13 - B.15. 
The effects of the inputs on a - AoA, appear in Figures B.18 - B.20. One should be 
noted to the scale of the damaged stabilator output - 10* in Figure 8.5 which is about 
three orders of magnitude smaller than the night stabilator response. In the damaged 
stabilator it is not surprising that the error is at the same magnitude as the actual output 
since the desired value is zero. The output of 10° is effectively zero. 

The capability of the network to react to individual input separately is even 
more emphasized in the damaged model. Although there is a totally different dynamics 
of the damaged input - left stabilator, as seen in Figure 8.5, it does not effect the training 
of the network to respond correctly to the rest of the inputs. Comparing Figures B.5- 


B.12 with B.13-B.20 one by one shows almost identical network behavior as was 
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expected from the similarities between the two models. Thus the observations on the 


undamaged aircraft hold to the damaged aircraft as well. 


Examining the results shown on Figures 8.4, 8.5 and Figures B.5-B.20 as a 


complete set describing the two networks used to emulate the undamaged and damaged 


aircraft separately, brings the following general observations. 


i 


The network resolves the relative influence of each input separately although 
while training in parallel the outputs included the cumulative contribution of all 
inputs. 


The network is able to resolve the similarities in the left and right control surface 
pairs, although each is treated separately. This is seen by the very close 
resemblance of the responses of each pair of inputs, which in the undamaged 
plant are identical, and also for the damaged plant except for the stabilators. 


The network is clearly able to identify the damage in the second plant as seen in 
Figures 8.5, B.13 and B.18 in comparison to Figures 8.4, B.5 and B.10. 


There is almost perfect agreement between the desired response and the network 


response. Practically they can be considered identical. For most cases the 
difference is less than 1%. 
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Figure 8.4: Normal Acceleration Response of the Undamaged Aircraft Due to 
Dipole Command in the Stabilators 
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Figure 8.5: Normal Acceleration Response of the Damaged Aircraft Due to Dipole 
Command in the Stabilators 
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B. EMULATION OF MULTIPLE PLANTS IN A SINGLE NETWORK 

For multiple-system emulation the same network structure as seen in Figure 8.1 is 
used for the results presented here, with p/J =1, p2=-1 for undamaged plant and p/ =-1, 
p2=1 for damaged plant. The network structure deliberately has had the same 
configuration as for a single plant emulation, so that comparisons can be made. Both the 
damaged and undamaged models were used to simultaneously train the single network 
via the procedure described in section IV.D. The results shown here are for the two 25" 
order models emulated together in a single network. It is important to mention that the 
procedure is general for any number of models, and two models were used here merely 


for demonstration, without any restrictions. 


1. Miultiple-Plant Emulation - Frequency Domain Results and discussion 

The examination starts with spectral analysis of the network compared to the 
desired model. Figure 8.6 shows the spectral response of the pitch rate q, the normal 
acceleration n, and the angle of attack AoA of the undamaged aircraft due to the right 
stabilator dst, (upper plots) and the left stabilator dst, (lower plots). Each graph shows 
the analytic frequency response in the solid line, which is the Bode plot of the outputs 
due to the corresponding input, together with the desired outputs in the dashed line, and 
network output values in the dash-dot line. Analogous plots for the leading edge and 
trailing edge control surfaces of the undamaged plant, appear in Appendix B, Figures 
B.21 and B.22 correspondingly. Figure 8.7 and the corresponding Figures B.23 and 
B.24 in Appendix B. presents the same information as Figures 8.6, B.21 and B.22, but 


for the damaged aircraft trained on the same network used for the undamaged aircraft. 
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The results from the spectral analysis obtained for the single-plant emulation 


for both undamaged and damaged aircraft hold for multiple plants emulated in a single 


network as well. The main results from Figures 8.6, 8.7, B.21 and B.22 show almost 


perfect match between the system Bode plot, the desired output and network output 


frequency response over the range of 0.03 Hz to 5 Hz. In Figure 8.7 the damage is 


clearly seen where there is no response of the aircraft to excitation due to the left 


stabilator (dstl) over the whole frequency range. 


Looking at Figures 8.6, 8.7 and Figures B.21-B.24 some general observations 


from the results are emerging. 


1. 


De 


The network learns the transfer function of each input individually. 


The unified network learns the full input-output transfer function of each plant 
individually. 


The network is able to resolve the similar input pairs such as left and nght 
stabilators, leading and trailing edge flaps, although each is treated separately. 
This is seen by the very close resemblance of the spectrums of each pair, nght 
and left inputs, which in the undamaged plant are identical, and also for the 
damaged plant except for the stabilators. This holds also for comparing the 
corresponding undamaged and damaged plots, i.e., the right stabilator part of 
Figures 8.6 and 8.7, as well as comparing Figure B.21 to B.23 and B.22 to B.24. 


The results of the multi-plant network are very similar compared to the results 
obtained from the single-plant network. Differences are due to the same 
reasoning given for the single-plant network in VII.A.1.a. 


The network is clearly able to identify the damage in the second plant as seen in 
Figure 8.7. 


There is a very good agreement between the network response and the analytical 
Bode plot over frequencies from 0.03 to 5 Hz. The differences are due to the 


insufficient persistent excitation in the input signal. 


The single network has learned the spectral behavior of the two plants as 
presented in Figures 8.6, 8.7, B.21-B.24. 
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Spectral Response of Undamaged Aircraft Due to Stabilators in a 
Multiple-Plant Single Network Emulation 
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Figure 8.7: Spectral Response of Damaged Aircraft Due to Stabilators in a 


Multiple-Plant Single Network Emulation 
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2. Miultiple-Plant Emulation - Time Domain Results and Discussion 

The time domain analysis shows the response of the outputs of the system, 
i.e., g, n, and a, to a dipole test signal driving each input. Due to the multiplicity of 
inputs and outputs, and in order not to disrupt the flow of reading, representative figures 
follow here, and the rest are gathered in Appendix B, and referred accordingly. The 
organization of all figures is identical to the ones shown and described here. 

Figurc 8.8 shows the normal acceleration response, n,, of the undamaged 
aircraft due to a dipole input command at the right stabilator dst, (upper plots) and at the 
left stabilator dst, (lower plots). The dotted line represents the dipole and is indicative 
of the switching times (the magnitude is +1 and is out of scale). The solid line is the 
desired behavior, and the dashed line in the neural net output. The response of the 
normal acceleration due to leading edge flaps dle. and dle,, and due to trailing edge flaps 
dte, and dte, appear in Appendix B, Figures B.28-B.29 respectively. The pitch rate 
response of the undamaged aircraft to a dipole command in all inputs is in Figures B.25- 
B.28. The corresponding angle-of-attack response appears in Figures B.30-B.32. 

The damaged aircraft was trained simultaneously on the same network. 
Figure 8.9 presents the normal acceleration response of the damaged aircraft due to the 
stabilators. The graphs of n, due to leading and trailing edge flaps appear in Figures 
B.36, B.37 respectively. The corresponding plots of the damaged aircraft pitch rate 
response to a dipole command in the various control surfaces are in Appendix B, Figures 


B.33 - B.35. The effects of the inputs on a@ - AoA, appear in Figures B.38 - B.40. One 
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should be noted to the scale of the damaged stabilator output - 10° in Figure 8.9 which 
is about three orders of magnitude smaller than the right stabilator response. 

Each graph shows the desired and network output for the command. The 
error between the two is presented in the lower plot of each pair. Note that the output 
errors are at the order of 10°. The relative error of input-output relations is from about 
5% down to 2%. It is interesting to see that the network which was trained by exciting 
all inputs simultaneously was equally trained on all inputs and both plants, regardless of 
the individual behavior and amplitudes. 

Examining the results shown on Figures 8.8, 8.9 and Figures B.25-B.40 as 
a complete set describing the single network used to emulate the undamaged and 
damaged aircraft together, brings the icin general observations. 


1. The results obtained for the single-plant network generally apply to the multiple- 
plant emulation as well as shown in Figures 8.8, 8.9, and Figure B.25-B.40. 


2. The output errors are greater than the errors achieved using a single plant 
emulation, and valued at about 2-5%. The differences are due to the greater 


complexity of the overall transfer function the network represents. 


3. The network is clearly able to distinguish the damaged plant from the undamaged 
model as of Figures 8.9, B.33 and B.38. 


4. There is almost perfect agreement between the desired response and the network 
response. Practically they can be considered identical. 


5. A single network is capable of emulating multiple systems with only a very small 
loss of accuracy due to the greater complexity. 
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Figure 8.8: Normal Acceleration Response of the Undamaged Aircraft Due to 
Dipole Command in the Stabilators in a Multiple-Plant Network 
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Figure 8.9: Normal Acceleration Response of the Damaged Aircraft Due to Dipole 
Command in the Stabilators in a Multiple-Plant Network 
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Figure 8.10 shows the normal acceleration response, n,, due to a random 


binary excitation of all inputs simultaneously. At the point indicated on the middle of 


the graph there is a transition from one model to the other. Observe that the time 


response is continuous in both cases. Each graph shows the desired and network output 


for the command five seconds before and after the switch. The error between the two 


is presented in the lower plot of each pair. Again plant #1 is the undamaged aircraft and 


plane #2 is the impaired model. The corresponding graphs for the pitch rate g and angle 


of attack a appear in Appendix B, Figures B.41 - B.42. 


The observations from the results in Figures 8.10, B.41 and B.42 are the 


following. 

1. There is a smooth continuous behavior of the output variables at the transition. 

2. The output error for a particular output is at the same level for both plants. 

3. The network is clearly able to distinguish the damaged plant from the undamaged 
model and respond accordingly. 

4. There is almost perfect agreement between the desired response and the network 
response when exciting all six control surfaces. This means that the network has 
learned the overall complex transfer function from the six inputs to the three 
Outputs altogether and not only per each individual input. 

5. The error at the transition point is slightly higher due to the regression, until the 


layer fills with data corresponding to the current model. Nevertheless, the error 
at the switch is at the same order of magnitude as in the normal operation. 


112 


nz Desired & NN Time response 














N ol. re NAY Deer 7 Se ee ene eine 7 
: : : 2 ! , ae | 
—_ — “tc ye 
- ! : 
6 7 8 9 10 11 | 72 ie) 14 eS) 
Time (sec] Do[-], NN[--] 
S 
us 
6 7. 8 9 10 14 eZ Le) 14 LS 
Time [sec] Do[=-], NN(--] @ 
1 nz Desired & NN Time response 
—™~ 
— 
16 Pr 18 19 20 24) 22 25 24 25 
Time [sec] DOo[-], NN(--] 
nz Error 
1,1 eel Wd 
M fat HP oe meena rag fh pay AEE op 
ia ial aT Tens a 
1 
: | lont #1: 
16 ez. 18 19 20 21 Pa. 25 24 pas 
Time [sec] DO[-], NN[--] af 


Figure 8.10: Normal Acceleration Response of the Single Network Emulating Two 
Plants Due to RBS Excitation of All Inputs 
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C. MINIMAL REALIZATION OF A SINGLE-NETWORK MULTIPLE-PLANT 

EMULATOR 

It has been shown that a single ANN is capable of emulating multiple-plants. A 
minimal realization of the multiple-plant network is now derived. The initial network 
shown in Figure 8.1 has two sub-networks, the forward on the left and the recurrent on 
the right. Each sub-network has two hidden layers with each hidden layer constructed 
of 36 PEs. This structure was successfully trained to emulate the two different aircraft 
models. This network served as the starting point for the minimal realization procedure. 

The regression layers, ComReg and OutReg of Figure 8.1, were already reduced 
based on observation of the decreasing valued of the regressed weights of more prior 
samples. The singular values of the weight matrices from the regression layers and 
through successive hidden layers were analyzed. Figure 8.11 shows the singular value 
plot of the initial weight matrix (36X36), connecting the two hidden layers of the 
recurrent sub-network. There are three dominant singular values and the rest are spread 
below 0.1. Close examination of Figure 8.11 shows that there are more significant 
singular values hiding in the initial range due to the random initialization. The number 
of PEs needed was estimated using equations (6.2), (6.3) with variation ratio of y=50Q0. 
The resultant layer needs only nine PEs. The singular values of the weight matrix 
connecting the same hidden layers in the minimal network are traced, as the network 
trains, in Figure 8.12. It can be seen that the same three dominant singular values 


appear in the minimal network, with the other six spread over the initial values. A lower 
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variation ratio would have predicted less PEs. It is the choice of the designer to set the 
threshold, with post-training examination of the achieved accuracy. 

A similar procedure was followed for the other three hidden layers. The singular 
values plots of the other weight matrices appear in Appendix C. On the recurrent side, 
the singular values of the weight between OutReg and outhl appear in Figures C.5 
before, and in Figure C.6 after training with similar behavior. 

In the forward sub-network, the singular values of the weight matrix between 
ComReg and comh! before minimization appear in Figure C.1, with values spread up to 
2.5. Using the same variation ratio, y=500, in equation (6.3) yields an estimation of 20 
PEs out of 36. The singular values of the minimized network are shown in Figure C.2, 
showing that the selection was quite cautious, and a smaller number might have been 
good enough for successful training. The Singular values of the weight matrix between 
comhI and comh2 are traced in Figures C.3 and C.4 respectively. 

The minimized network appear in Figure 8.13. The squares represent the 
processing elements. According to the current input value and the weights connecting 
to a PE the output value is viewed by size and color, converted here to shades of grey. 
It is seen that the activation from one layer to another is transferred differently by the 
PEs. This is clearly seen in ComReg and comh! layers. Given an input in terms of 
control commands at time (t=k) at the Jn layer of Figure 8.13, and a model selection at 
the input nodes denoted p/ or p2, the network is trained to generate an output of the 


selected aircraft model at time (¢=k+/). 
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Looking at the singular values plots in Figures 8.11, 8.12 and C.1-C.6, the 


following observations are drawn. 


le 


Singuiar Values 


The singular values plot express the number of PEs carrying the activation levels 
between two successive layers. 


The singular values plots trace the dynamics of the training, and may supply an 
estimation on when the network is trained, by the fact that the singular values 
converged to their final values. The slopes of the curves need to be divided by 
the temporal learning rate to show convergence. Figure 8.11 reflect convergence 
after about 500k iterations, whereas in the minimized network, as seen in Figure 
8.12, after about 800k the network is not well converged yet. 


The plots should be carefully examined for small but important singular values 
that might hide in the initial randomization values. This has happened i.e., in 
Figure C.3 and extracted also using equations (6.2) and (6.3). 


The network should be trained long enough for the number of significant singular 
values to become evident, depending of course on the learning rate and initial 
randomization. This phenomenon is clearly seen in Figure 8.11 and 8.12 where 
the third dominant singular value emerge after about 100k iterations. 


The singular values plots act as a feature detector, and can give physical 
explanation to the training characteristics of the network. In the current example 
the third dominant singular value in Figures 8.11, 8.12 C.5 and C.6, probably 
can be associated with the phugoid mode, which needed more iterations to have 
a significant number of presentations (about 10) for the network to respond to that 
mode. Once the network sensed this mode the adaptation of the weight to include 
the phugoid mode information was rather quick. The short period mode was 
probably related to the quick rising singular values. 
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Figure 8.11: Singular Values Plot of the Non-Minimized Network 
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Figure 8.12: Singular Values Plot of the Minimized Network 
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Figure 8.13: Minimal Realization of a Network for Multiple-Plant F/A-18 Emulation 
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1. Minimal Realization of a Multiple-Plant Network - Frequency Domain 

Results and Discussion 

The spectral analysis compares the network with the desired model and with 
the performance of the non-minimized network. Figure 8.14 shows the spectral response 
of the pitch rate g, the normal acceleration n, and the angle of attack AoA of the 
undamaged aircraft due to the night stabilator dst, (upper plots) and the left stabilator dst, 
(lower plots). Each graph shows the analytic frequency response in the solid line, which 
is the Bode plot of the outputs due to the corresponding input, together with the desired 
outputs in the dashed line, and network output values in the dash-dot line. Analogous 
plots for the leading edge and trailing edge control surfaces effects on the undamaged 
plant, appear in Appendix B, Figures B.43 and B.44 correspondingly. Figure 8.15 and 
the corresponding Figures B.45 and B.46 in Appendix B, presents the same information 
as Figures 8.14, B.43 and B.44, but for the damaged aircraft trained on the same 
minimized network used for the undamaged aircraft. 

The results from the spectral analysis obtained for the multiple-plant 
emulation of both undamaged and damaged aircraft in a single network, summarized in 
section B, hold for the minimized network emulator as well. The main results from 
Figures 8.14, 8.15, B.43 and B.44 show almost perfect match between the system Bode 
plot, the desired output and network output frequency response over the range of 0.03 
Hz to 5 Hz. In Figure 8.15 the damage is clearly seen where there ts no response of the 
aircraft to excitation due to the left stabilator (dstl) over the whole frequency range. The 


fluctuations seen in Figure 8.15 are due to initial randomization, numerical noise and the 
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random nature of the test signal. The values are at the order of -80dB, which can be 
neglected. 

The results here are to be compared to the multiple-plant emulation in Figures 
8.6,8.7 and Figures B.21 - B.24 correspondingly. The combined network is representing 
now a minimal network realization of the non-linear functional relationship between the 
inputs and outputs of the network that describes the damaged and undamaged aircraft. 

Looking at Figures 8.14, 8.15 and Figures B.43-B.44, some general 
observations from the results are drawn. 

1. The observations for the Multiple-plant emulation generally apply to the minimal 
realized multiple-plant emulation as well. There are instances that the minimized 
network shows even better results over the spectrum, like in the leading edge 
surfaces effect on both undamaged and damaged aircraft, as seen in comparing 
Figures B.21 for the multiple-plant emulation, and in Figure B.43 for the minimal 
realization. The improved behavior is probable due to the elimination of excess 


PEs and excessive connections in the minimized neural network. 


2. The single minimized network has learned the spectral behavior of the two plants. 
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Figure 8.14: Spectral Response of Undamaged Aircraft Due to Stabilators in a 
Minimally Realized Multiple-Plant Single Network Emulator 


120 


Fjont #2 Freq. response of q due to dstr 








Gain dB 
Gain dB 


Bb ; aig eee 
10-3 10-2 10m 109 





10! 10-5 10-2 LO=" 10° 10! 


Frequency (Hz) Frequency (Hz) 
Pignt #2 Freq. response of AoA due to dstr 10 Plant #2 -— LEGENO 








Desired Output 
m -—20 
cel NN Qutput 
— 
S -40 Sys Bode 
a ae 
10-5 10-2 1e-" 19° 10! 0 5 10 
Frequency (Hz) ay 
Bont #2 Freq. response of q due to dst Plont #2 Freq. 
eal eet es) 
ie On vs ee, 
5 fs] 
a) ae ° 
-100 MA were =106 7 
10-4 1072 es! 10° 10! 10-5 10-2 10>! 199 ee 
Frequency (Hz) Frequency (Hz) 
Pignt #2 Freq response of AoA due to dst 10 Plant #2 — LEGENO 


4B Desired Output 

fee) — 

. NN Qutput 

& —100 Sys Bede 
—150 


Oz? TOs TO sa 10° 10) 


Frequency (Hz) 





3mf 


Figure 8.15: Spectral Response of Damaged Aircraft Due to Stabilators in a 


Minimally Realized Multiple-Plant Single Network Emulator 
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2. Minimal Realization of a Multiple-Plant Network - Time Domain Results 

and Discussion 

The time domain analysis compares the response of the outputs of the system, 
.1e., g, m, and a, with the network reaction to a dipole test signal driving each input. 
Figure 8.16 shows the normal acceleration response, n., of the undamaged aircraft due 
to a dipole input command at the right stabilator dst, (upper plots) and at the left 
stabilator dst, (lower plots). The dotted line represents the dipole and is indica.ive of the 
switching times (the magnitude is +1 and is out of scale). The solid line is the desired 
behavior, and the dashed line in the neural net output. The response of the normal 
acceleration due to leading edge flaps dle, and dle,, and due to trailing edge flaps dte, and 
dte, appear in Appendix B. Figures B.50-B.51 respectively. The pitch rate response of 
the undamaged aircraft to a dipole command in all inputs is in Figures B.47-B.49. The 
corresponding angle-of-attack response appears in Figures B.52-B.54. 

The damaged aircraft was trained simultaneously on the same minimized 
network. Figure 8.17 presents the normal acceleration response of the damaged aircraft 
due to the stabilators. The graphs of n, due to leading and trailing edge flaps appear in 
Figures B.58, B.59 respectively. The corresponding plots of the minimized network 
representation of the damaged aircraft pitch rate response to a dipole command in the 
various control surfaces are in Appendix B, Figures B.55 - B.57. The effects of the 
inputs on a@ - AoA, appear in Figures B.60 - B.62. One should be noted to the scale of 
the damaged stabilator output - 10° in Figure 8.17 which is about three orders of 


magnitude smaller than the right stabilator response. 
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Each graph shows the desired and network output for the command. The 
error between the two is presented in the lower plot of each pair. Note that the output 
errors are at the order of 10°. The relative error of input-output relations is about 5% 
at most down to 1%. The network was trained by exciting all inputs simultaneously and 
equally performs on all inputs of both plants, regardless of the individual behavior and 
amplitudes. 

Examining the results shown on Figures 8.16, 8.17 and Figures B.47-B.62 
as a complete set describing the single minimal network used to emulate the undamaged 
and damaged aircraft together, and comparing to the non-minimal single network for 
multiple-plant emulation in Figures 8.8,8.9 and Figures B.25-B.40 correspondingly, 
brings the following general observations. | 


1. The results obtained for the single network of multiple-plant generally apply to 
the minimally realized single network of multiple-plant emulation as well. 


2. The output errors are at the same level as the errors achieved using a mom- 
minimal single network for multiple-plant emulation, and valued at about 1-5%. 
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Figure 8.16: Normal Acceleration Response of the Undamaged Aircraft Due to 
Dipole Stabilators Command in a Minimized Multiple-Plant Network 
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Figure 8.17: Normal Acceleration Response of the Damaged Aircraft Due to Dipole 
Stabilators Command in a Minimized Multiple-Plant Network 
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Figure 8.18 shows the normal acceleration response, n., due to a random 
binary excitation of all inputs simultaneously. At the point indicated on the middle of 
the graph there is a transition from one model to the other. Each graph shows the 
desired and network output for the command five seconds before and after the switch. 
The error between the two is presented in the lower plot of each pair. Plant #1 is the 
undamaged aircraft and plant #2 is the impaired model. The corresponding graphs for 
the pitch rate q and angle of attack a appear in Appendix B, Figures B.63 - B.64. The 
results here are to be compared to the non-minimal single network for multiple-plant 
emulation in Figure 8.10 and Figures B.41 - B.42 correspondingly. 

The observations from the results in Figures 8-18, B.63 and B.64 in 
comparison to the results of the non-minimal network in Figures 8.10 and B.41, B.42 
are the following. 


1. The observations for the non-minimal single network for multiple-plant emulation 
also hold for the minimal realization case. 


2. There is no degradation in performance due to the less complex network achieved 
by the minimal realization procedure. 
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Figure 8.18: Normal Acceleration Response of a Minimally Realized Single Network 
Emulating Two Plants Due to RBS Excitation of Al! Inputs 
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D. EMULATION ANALYSIS 
The capabilities of an artificial neural network to emulate MIMO dynamical 
systems have been shown. The main results are the following. 


1. A combination of feedforward and recurrent networks is used to emulate MIMO 
dynamical systems. 


2. The network structure is selected to match the class that the system belongs to. 


3. A single network is capable of emulating multiple plants under proper training 
scheme. 


4. More presentations are needed to emulate multiple plants. 
5. A minimal realization of a network can be found, with no degradation in 


performance. Corollary: there is no unique solution to the emulation / 
identification problem. 
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E. ADAPTIVE CONTROL OF MIMO DYNAMICAL SYSTEM 

The control results include two parts. Initially a model follower controller is 
designed for each of the undamaged and damaged aircraft separately. Then the two 
controllers are combined into a single network. The unified network is an adaptive 
controller. Based on the input information, the network generates the required control 
actions such that the plant output will behave similar to the reference model. The input 
information includes the pilot command p,, the plant output feedback, and prior 
measurements of plant input (network output) and output as indicators of plant type, of 


either undamaged or damaged aircraft. 


1. Model-Reference Controller - Frequency Domain Results and Discussion 

The spectral response of the pitch rate g, the normal acceleration n, and the 

angle of attack AoA of the undamaged (plant #1) and damaged (plant #2) aircraft due to 
pilot commands p, is given in Figure 8.19. Each graph shows the analytic frequency 
response - the solid line, which is the Bode plot of the reference system outputs due to 
the input, together with the desired values and network output values, both analyzed 
using FFT, while exciting the input (pilot) command with random binary sequence 
(RBS). The desired value shown by the dashed line is calculated using the reference 
model, and the dash-dot line is the neural network output. Figure 8.20 presents the 
same information using the unified controller, which is capable of controlling both the 
undamaged and damaged plants simultaneously. The normalization of the models to +1 


cause the shift of about -15dB from the zero dB line. 
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The observations from the spectral results in Figures 8.19 and 8.20 are 


the following. 


le 


Each plant is controlled such that the overall system behave very close to the 
reference model over the entire spectral range of interest. The difference at the 
low frequencies is a problem of persistent excitation. For piloted aircraft control 
usually the indicated difference at the low end does not present a control problem. 
The human pilot compensate for the low frequency dynamics with no accessional 
control effort. 


The spectral behavior of the unified controller in Figure 8.20 is similar for both 
models, with good agreement to the analytical reference, although less accurate 
than the single controller shown in Figure 8.19. The difference between the 
unified and single controllers is due to the greater complexity of the unified 
controller. 
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Figure 8.19: Spectral Response of Undamaged and Damaged Controlled Aircraft 
Due to Pilot Longitudinal Stick Command, Single Controller 
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Figure 8.20: Spectral Response of Undamaged and Damaged Controlled Aircraft 
Due to Pilot Longitudinal Stick Command, Unified Controller 
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2. Model-Reference Controller - Time Domain Results and Discussion 

The time domain analysis shows the response of the outputs of the system to 
a dipole test signal at the pilot command p,. The dipole presents a one second full stick 
movement that demand a 5g acceleration from the aircraft, in both directions. Each 
graph shows the reference, which is the desired output, and the network output q, n, and 
a for the dipole command. The error between the two is presented in the lower plot of 
each pair. Figure 8.22 presents the same information on the damaged aircraft, each 
controlled with a separate neural network controller. The stabilators, leading and trailing 
edge flaps actuators commands in response to a dipole input are shown in Figure 8.23 
for the undamaged aircraft, and in Figure 8.24 for the damaged aircraft. 

The performance of the unified controller in the time domain is tested by 
presenting a full stick RBS input as the pilot command. The response of the network 
compared to the reference model is presented in Figure 8.25 with the corresponding 
actuators commands in Figure 8.26. 

Each graph shows the desired and network output for the command. The 
error between the two is presented in the lower plot of each pair. Note that the output 
errors are at the order of 10° or less. The relative error of all input-output relations is 


about 1-5%. 
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Examining the results shown on Figures 8.22 - 8.26 as a complete set 


describing the neural network controlled undamaged and damaged aircraft, brings the 


following general observations. 


1. 


A neural network based controller was designed for the undamaged as well as the 
damaged aircraft. 


The network finds a suitable solution, although not unique. In the undamaged 
aircraft the three control pairs are symmetric. Depending on the initial weights 
randomization the network is trained to find an equilibrium solution which is not 
necessarily symmetric as one could expect in Figure 8.23. 


The network is clearly able to identify the damage in the second plant and 
disregard the malfunctioning surface (the left stabilator command is practically 
zero). The neural network learns to compensate for the loss of control surface 
with the rest of the actuators, mainly with the adjacent right stabilator as can be 
clearly seen in Figures 8.24 and 8.26. 


The unified controller is able to identify the change in the plant and switch over 


to generate the required control commands. The physical output variables remain 
continuous, while the behavior of the control surfaces changes dramatically. 


Control Analysis 


In this section the capabilities of an artificial neural network to adaptively 


control a time-varying MIMO dynamical systems have been shown. The main results 


are the following. 


I. 


An artificial neural network was trained as controller to drive a plant in a model 
following scheme. 


The neural network controller was able to compensate for a loss of control action. 


A single network is capable of controlling multiple plants in an adaptive fashion. 
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Figure 8.21: Time Response of Undamaged Controlled Aircraft Due to Pilot Dipole 
Stick Command, Single Controller 
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Figure 8.22: Time Response of Damaged Controlled Aircraft Due to Pilot Dipole 
Stick Command, Single Controller 
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Figure 8.23: Actuation Commands of Undamaged Controlled Aircraft Due to Pilot 
Dipole Stick Command, Single Controller 
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Figure 8.24: Actuation Commands of Damaged Controlled Aircraft Due to Pilot 
Dipole Stick Command, Single Controller 
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Figure 8.25: Time response of Controlled Time Varying Aircraft Due to Pilot RBS 
Stick Excitation, Unified Controller 
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Figure 8.26: Actuation Commands of Controlled Time Varying Aircraft Due to 


Pilot RBS Stick Excitation, Unified Controller 
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IX. CONCLUSIONS AND RECOMMENDATIONS 


This work considered dynamical systems in association with artificial neural 
networks. The basic theory of neural network was presented. Methods for emulation 
and control of dynamical systems using neural networks were developed. The important 
results are summarized here and some recommendations for further research are 


highlighted. 


A. MAIN RESULTS 

This research established a general approach for identification and emulation of 
non-linear MIMO time-varying dynamical system. Four general MIMO dynamical 
models are presented. Based on the selected model a corresponding network structure 
was introduced. The structure associated a combination of feedforward and recurrent 
artificial neural networks. A serial-parallel training scheme was used to train the 
network. A single network was shown to be capable of representing multiple dynamical 
systems from the same class, under the balanced random training procedure. The model 
switching in the procedure did not affect the accuracy of the network. 

An algorithm to find minimal realization of a neural network was also presented. 
It was shown that there was no loss in performance of the reduced neural network, 
comparing to a more complex network that does the same task. 

Further in the research a neural network was trained to control a MIMO dynamical 


system using a back-propagation approach through the plant training procedure, which 
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is a variation of ite basic back-propagation algorithm. The network was trained in a 
model-follower control architecture without a direct access to the desired values and 
without the need for specific knowledge on the controlled system. A single network 
controller was trained to represent multiple controllers, thus forming an adaptive 
controller. No estimation process was involved in the procedure. 

Artificial neural networks were found to be a powerful tool for system 


identificatior. and in designing adaptive controllers. 


B. RECOMMENDATIONS FOR FURTHER RESEARCH 

Several topics that were addressed in the research need to be further investigated. 
The issues are discussed below with possible directions of approach. 

In the minimal realization procedure, singular values of a weight matrix were 
involved. In the LTI system emulation the singular values had physical meaning. It 
seems that the singular values may give more information about the system. It can also 
serve as an indicator of network training progress and also as a feature detector. The 
singular value analysis of the weight matrix needs more investigation in the above 
directions as well as in improving the minimization procedure to more definite results. 

Training recurrent networks has still many open aspects. Since recurrent networks 
are vital in association with dynamical system, more research is needed in the theoretical 
and practical aspects of the training. Methods for stable algorithms of training hidden 
recurrent layers are affiliated with this direction of research as well. General constraints 


on the learning rate for stable training are greatly desired. Proper learning rate may also 
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improve convergence rate. Methods involving numerical analysis in stability theory may 


be an approach. 
Algorithms for utilizing adaptive neural networks controller in real time and on-line 


applications may be seen as a direct extension of this work. This area of research may 


fined many current applications, especially in the non-linear regime. 
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APPENDIX A. F/A-18A MATHEMATICAL MODEL 


The F/A-18A mathematical model for the damaged and undamaged configurations 
is given this MATLAB™ script file. The aerodynamical stability derivatives are shown 


in the body of the program. 


*% f18longN.m 

% 

% This file calculates the continuous state Space representation of the 
% longitudinal modes of the f18 (with or without malfunctions) 

% fighter aircraft. 

% Then the matrices are balanced and discretized according to a given 
% sampling interval - ts. The results are to be loaded to file: 

% f18long_model .h 

% The state vector: x=[(u(ft/sec),w(ft/sec),q(rad/sec), theta(rad)) ’ 
% The control vector: del=[(dstr,dstl,dler,dlel,dter,dtel]’ (deg) 

% The calculation assumes that the gusts are negligible. 

% Trim AOA=alpha0=arctan(W/U), and theta0=gamma0+alpha0. 

% This uses cscaleb.m 

cs=[); 

adirectory='/f11669/home/dror/f18/Emul/’; 
directory='/f11669/home/dror/f1i8/Cntrl/’; 


ts=0.02; % (sec)- Sampling interval 

al pha0=2.6184*pi/180; % [rad] - Trim angle of attack 
gamma0=0.0; % [rad] - Path angle 
theta0=alpha0+gamma0; % (rad) - a/c body angle 
U=646.42; % (ft/sec) - True airspeed 
alt=10000; % (ft) - altitude 

M=0 .6; % (#] - Mach number 

g=32.174; % (ft/sec*2]) - gravity constant 
W=U* tan(alphaQ); 

W=0; % In trim axes 

“theta0=gamma0; % In trim axes 


disp(’ Definitions of dimensional stability derivatives’ ) 


% Name Value Units Partial deriv. of With respect to 
yA 

XU= -0.13257e-1; 4(1/sec] Longitudinal force Forward velocity 
XwW= 0.71265e-1; a(1/sec] Longitudinal force Vertical velocity 
xq= 0.0: “”l(ft/(rad-sec)] Longitudinal force Pitch rate 

xwd= 0-0. A) Longitudinal force Vertical acceler 
xdsb= 0.0; 4(ft/sec*2)/rad Longitudinal force Speed brake 

xdth= 0.14257; 4 (pcet/sec*2) Longitudinal force Throttle 

xds= 0-0; ’(ft/sec*2)/rad Longitudinal force Horiz stabilator 
xdlf= 0.0; a(ft/sec*2)/rad Longitudinal force Lead edge flap 
xdtf= 0.0: *(ft/sec*2)/rad Longitudinal force Trail edge flap 
Zu= “0. /5957e- 1: 4(1/sec] Vertical force Forward velocity 
Zw= =enloeo 4(1/sec] Vertical force Vertical velocity 
zq= 5) Aes Ip) 5 a(ft/C(rad-sec)] Vertical force Pitch rate 

zwd= -0.006917; ’(--] Vertical force Vertical acceler 
zdsb= 0.0; “a(ft/sec*2)/rad Vertical force Speed brake 

zdth= 0.0; 4(pet/sec*2] Vertical force Throttle 

zds= -129.5; a(ft/sec*2)/rad Vertical force Horiz stabilator 
zdif= 19.43; “(ft/sec*2)/rad Vertical force Lead edge flap 


144 


zdtf= -149.0; a(ft/sec*2)/rad Vertical force Trail edge flap 


mu= -0.12988e-4; 4(1/sec- ft] Pitching moment Forward velocity 

mw= =Ooilsste-1; 4(1/sec- ft] Pitching moment Vertical velocity 
= -0.59346; 4(1/sec] Pitching moment Pitch rate 

mwd= -0.34049e-3; “C1/Ft) Pitching moment Vertical acceler 

mdsb= 0.0; 4(1/sec*2)/rad Pitching moment Speed brake 

mdth= 0.0; “(pcet/(sec*2-ft)JPitching moment Throttle 

mds= -15.6; %(1/sec*2)/rad Pitching moment Horiz stabilator 

mdlif= -1.609; 4(1/sec*2)/rad Pitching moment Lead edge flap 

mdtf= 1.499; 4(1/sec*2)/rad Pitching moment Trail edge flap 

% i ee ee ee 


% Actuators are scaled to +/-1 


% STABILATOR 
stln=(1/(82.9)*2, 2*(0.068)/(82.9), 11;% left side numerator 
st ld=conv( (1/(36.4)*2,2*(0.41)/(36.4), 1), (1/¢105 .3)*2, 2*(0.59)/105.3,1]); 


mxstl=10.5; % maximum positive deflection of left stabilator [deg] 
mnstl=-24; % minimum negative deflection of left stabilator [deg] 
strn=stln; % Tight side numerator 

strd=stld; *% right side denominator 

mxstr=10.5; % Maximum positive deflection of right stabilator [deg] 
mnstr=-24; 4 minimum negative deflection of right stabilator [deg] 


if(model=='2') stln=0; end; 
(Fstl,Gstl,Hstl,Ostlj=tf2ss(stin,stld); 
(Fstr,Gstr,Hstr,Dstrj=tf2ss(strn,strd); 


nust l=(mxstl-mnstl)/2; 4 Left stabilator scaler 
bist l=(mxstl+mnstl)/(2*nustl); % Left stabilator bias 
nustr=(mxstr-mnstr)/2; 4% Right stabilator scaler 
bistr=(mxstr+mnstr)/(2*nustl); % Right stabilator bias 


% LEADING EDGE FLAP 


leln=1; 

leld=conv( (1/(26.9), 1],(£1/¢82.9), 1]); 

mxlel=34; % Maximum positive deflection - left leading edge flap [deg] 
mnlel=-3; % minimum negative deflection - left leading edge flap [deg] 


lern=leln; 

lerd=leld; 

mxler=34; % Maximum positive deflection of right leading edge flap [deg] 
mnler=-3; 4% minimum negative deflection of right leading edge flap [deg] 
(Flel,Glel,Hlel,Dlelj=tf2ss(leln, leld); 

(Fler,Gler,Hler,Dlerj=tf2ss(lern, lerd); 


nulel=(mxlel-mnlel)/2; % Left leading edge flap scaler 
bilel=(mxlel+mnlel)/(2*nulel); % Left leading edge flap bias 
nuler=(mxler-mnler)/2; % Right leading edge flap scaler 


biler=(mxler+mnler)/(2*nulel): % Right leading edge flap bias 


% TRAILING EDGE FLAP 


teln=1; 

teld=(1/(35)*2, 2*(0.71)/(35), 1]; 

mxtel=45; % Maximum positive deflection -left trailing edge flap [deg] 
mntel=-8; % minimum negative deflection -left trailing edge flap [deg] 
tern=teln; 

terd=teld; ° 

mxter=45; % maximum positive deflection -right trailing edge flap(deg] 
mnter=-8:; 4% minimum negative deflection -right trailing edge flap[deg] 


(Ftel,Gtel,Htel ,DteljJ=tf2ss(teln,teld); 
(Fter,Gter,Hter ,Oterj=tfess(tern,terd); 


nutel=(mxtel-mntel)/2>; 4% Left trailing edge flap scaler 
bitel=(mxtel+mntel )/(2*nutel ); 4 Left trailing edge flap bias 
nuter=(mxter-mnter )/2; % Right trailing edge flap scaler 
biter=(mxter+mnter )/(2*nutel); % Right trailing edge flap bias 
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% Scaling and biasing 

NU=(nustl nustr nulel nuler nutel nuter]; % scale input from+/-1 to physical 
bi=(bistl bistr bilel biler bitel biter]; % input shift (bias) 

NU=max(abs( ((mxstl;mnstl), (mxstr;mnstr), (mxlelsmnlel), (mxler;mnler]), 

(mxtel smntel), (mxter;mnter]])); 

bi=zeros(bi); 

disp(’ sensors transfer functions, and their state space representation’ ) 

yA +—+}_}_}_$—}-_+- $+} _}-_}-_$ $45} _$- $$} _}-_}—}—-}— $f — ff et tt ht lt — 


% RATE GYRO 

raan=(1/C131.7) 1); % Pitch rate rate-gyro 
rgqd=conv( (1/( 167.8), 1),(1/¢€461.7), 1]); 

(Frgq, Grgq, Hrgq,Drgq) =tf2ss(rgqn, rgqd); 


% ACCELEROMETER 

aczn=(1/(235.8), 1]; 4% Normal accelerometer 
aczd=(1/(395.3)*2, 2*(0.96)/(395.3), 1); 

(Facz,Gacz,Hacz,Dacz] =tf2ss(aczn, aczd); 


% ANGLE OF ATTACK 

aan=1; 

aad=(1/(14), 1); 

(Faa, Gaa, Haa,Daa) =tf2ss(aan, aad); 


disp(’ Calculation of longitudinal matrices’ ) 
% calculating continuous longitudinal dynamic matrix. 
% In:Xx=(u,w,q,theta]’, Out:Xxd 


Fx(2,:)=(zu zw (zqtU) -g*sin(theta0)]/(1-zwd); 
Fx(1,:)=(xu xw (xq-W) -g*cos(theta0))] +xwd*Fx(2,:); 
Fx(3,:)=(mu mw mq 0 ] +mwd*Fx(2,:); 
Fx(4,:)=(0, 0, 1, 0 a 


% calculating continuous longitudinal control matrix. 
% In:Ox=(dstx,dlex,dtex]’, Out:Xxd 


Gx(2,:)=(zds zdlf zdtf]/(1-zwd); 
Gx(1,:)=(xds xdlf xdtf]+Gx (2, :)*xwd; 
Gx(3,:)=(mds mdl f mdtf]+Gx(2, :)*mwd; 
Gx(4,:)=(0, 0, O J]; 


% calculating longitudinal measurement matrix. In:Xx, Out:Yx=(q,nz, alpha] ’ 
Hx(1,:)=(0, 0, 1, 0); 


Hx(2,:2)=Fx(2,:)-Hx(1,:)*U; %NZ=wdot -U*q=U*( al phadot-q) 
Hx(3,:)=(0, 1/U, 0, OJ; 


% calculating longitudinal control feedforward matrix. 
% In:Dx=(dstx,dlex,dtex)’, Out:Yx 
Oxcl..* }=(0-20. O]- 


Ox(2,:)=Gx(2,:); 
Ox(3,:)=(0, 0, 0); 


% calculating longitudinal controls distribution. 
4 In:0=(dstr,dstl,dler,dlel,dter,dtel) ’ ,Out:Dx=(dstx,dlex,dtex] ’ 


lng=zeros(3,6); 
lng(1,1)=0.5; lng¢(1,2)=0.5; lng(2,3)=0.5; lng(2,4)=0.5; lng(3,5)=0.5; lng(3,6)=0.5; 
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disp(’ appending system matrices’ ) 

% system dynamic matrix. In:X=[u,w,qg,theta]’, Out: Xd. 

% control matrix. In:O=[dstr,dstl,dler,dlel,dter,dtel]’, Out:Xd 
% measurement matrix. I[n:X, Out:Y=[g,nz,alpha]’. 

% control feed-forward matrix. In:0, Out:Y . 


Fm=F x; 
Gm=Gx* lng; 
Hm=Hx ; 
Om=0x* lng; 


% convert from (per radian) to (per degree) and from (ft/sec*2) to (g) 
disp(’convert (/rad) to (/deg) & (ft/sec*2) to (g)’); 

[nx, nu] =size(Gm); [ny, nu] =size(0m); 

r2d=180/pi; 

d2r=1/r2d; 

Nu=eye(nu)*dér; % Get inputs in degrees 
Nx=diag((1,U,1,1]); % change w to alpha 

Ny=diag( (d2r,g,d2r]); % convert to degrees and g’s 
Fm=Nx \Fm*Nx > 

Gm=Nx\Gm* Nu; 

Hm=Ny\Hm*N x 3 

Om=Ny\0m*Nu; 

end; 


disp(’ actuators state variable model’) 


([FA,GA,HA,DA]=append(Fstr,Gstr,Hstr,Ostr,Fstl,Gstl,Hstl,Ostl); 
[FA,GA,HA,DA] =append(FA,GA,HA,OA,Fler,Gler,Hler,Oler); 
[FA,GA,HA, 0A] =append(FA,GA,HA,DA,Flel,Glel,Hlel,Dlel); 
[FA,GA,HA,OA]=append(FA,GA,HA,DA, Fter,Gter,Hter,Oter); 

[FA, GA, HA,DA] =append(FA,GA,HA,DA,Ftel,Gtel,Htel,Otel); 

% Scale inputs 

GA=GA*diag(NU); DA=DA*diag(NU); 


disp(’ sensors state variable model’) 


(FS,GS,HS,0S] =append( Frgq, Grgq, Hrgq,Orgaq, Facz,Gacz,Hacz,0acz); 
(FS,GS,HS,0S] =append(FS,GS,HS,0S,Faa,Gaa,Haa,Daa); 


disp(‘check and balance matrices’ ) 

(Fm,Gm,Hm]=ckbal(‘ Fm’ ,’Gm’,’Hm! ,Fm,Gm, Hm); 
[FA,GA,HA]=ckbal(’FA’, ‘GA’, ‘HA’, FA,GA,HA); 
[FS,GS,HSJ=ckbal(’FS’,’GS’,’HS’,FS,GS,HS); 


disp( ‘cascading the system’ ) 
([Fps,Gps,Hps ,Ops]=series(FA,GA,HA,OA, Fm, Gm, Hm,0Om); 
([Fps,Gps,Hps ,Ops]=series(Fps,Gps,Hps,Ops,FS,GS,HS,0S); 


% must discretize AFTER cascading continuous systems 
disp(‘discretize the matrices’) 
[Fod, God] =c2d(Fps,Gps,ts); 


4% find the extrimum values of each output due to each controller and 
% normalize to +/- 1. 

by=zeros(1,ny); 

cs=ones(1,ny); 

11=input([/D0 you wish to normalize outputs to +/-1? [y]l/n ‘],'s‘); 
if isempty (ii), 11=’y’ ;end; 

if ij==s/y! “Normalize 
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disp(‘Normalizing...’) 
format compact 
if(model~='1') eval( (‘load ‘ directory ‘zcsi.mat‘]);% get cs,by 
cscaleb(‘d’ , Fpd,Gpd,diag(cs)\Hps,diag(cs)\Dps,bi,by,ts,1,1); 
by,cs % show by and cs 
else 
l=input((’Enter # of iterations (le6]): ‘]); 
if isempty(l), l=1e6;end; 
while 1 
(cs, by]=cscaleb('d’ ,Fpd,Gpd,Hps,0ps,bi,by,ts,1,1); 
i=input(/Do you want to recalculate cs, by? y/{(n) ‘,’s‘); 
if (isempty(i)) break;end; . 
if(iz='n’) break; end; 
end; 
izinput(’Do you want to edit cs, by (via maxy,miny)? y/(n] ‘,’s’); 
if Cisempty(i)) 1=!’n’;end; 
if(iz='n’) ; 


else 
keyboard; 
by=-0.95*(maxy+miny) ./(maxy-miny) 
cs= (maxy-miny)/(2*0.95) % The 0.95 is 5% safety factor 
end; 
format; 
end; 
Hps=diag(cs)\Hps; 
Ops=diag(cs)\Dps; 
end; 
i, —Trrrerrrrerrererrerrrrrrrrrrrrrere eee ere eee eee 


i=input((/Do you want to save matrices to ’ directory ‘ for NN use? [y])/n ‘J],‘s‘); 
if isempty(i), 1=’y’;end; 

if i==’y’% Save 

eval(({’save ‘ directory ‘zFld’ model ’.mat Fpd’)); 

eval( (‘save ‘ directory ‘zGld’ model ‘’.mat Gpd’)); 

eval( (‘save ’ directory ‘ZHld’ model ‘.mat Hps’)); 

eval( (‘save ' directory ‘zDld’ model ‘’.mat Ops’)); 

eval( (‘save ’ directory ‘zBil’ model ‘.mat bi’)); 

eval(({’save ‘ directory ‘zBol’ model ‘.mat by’])>; 

eval((’save ‘ directory ‘zcs’ model ‘.mat cs by NU Nu Nx Ny’)); 
end; 

format 
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APPENDIX B. TIME AND FREQUENCY RESPONSES OF 
ARTIFICIAL NEURAL NETWORK EMULATION OF 
DAMAGED AND UNDAMAGED F/A-18A AIRCRAFT 


The graphs in Appendix B are referenced in chapter VU. The reader is 
referred to chapter VIL for comprehensive explanation and related discussion 


regarding the plots that appear in this appendix. 
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Spectral Response of Undamaged Aircraft Due to Leading Edge Flaps, 
Single Plant, Single Network Emulation 
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Figure B.2: Spectral Response of Undamaged Aircraft Due to Trailing Edge Flaps, 


Single Plant, Single Network Emulation 
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Spectral Response of Damaged Aircraft Due to Leading Edge Flaps, 
Single Plant, Single Network Emulation 
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Figure B.4: Spectral Response of Damaged Aircraft Due to Trailing Edge Flaps, 


Single Plant, Single Network Emulation 
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Figure B.5: Pitch Rate Response of Undamaged Aircraft Due to Stabilators Dipole 
Command, Single Plant, Single Network Emulation 
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q Desired & NN Time response due to dier Plant #1 
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Figure B.6: Pitch Rate Response of Undamaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Single Plant, Single Network Emulation 
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q Desired & NN Time response due to dter Plant #1 
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Figure B.7: Pitch Rate Response of Undamaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Single Plant, Single Network Emulation 
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nz Desired & NN Time response due to dier Plant #1 
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Figure B.8: Normal Acceleration Response of Undamaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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nz Desired & NN Time response due to dter Plont #1 
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Figure B.9: Normal Acceleration Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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AoA Desired & NN Time response due to dstr Plant #1 
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Figure B.10: Angle of Attack Response of Undamaged Aircraft Due to Stabilators, 
Dipole Command, Single Plant, Single Network Emulation 
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AoA Desired & NN Time response due to dier Plant #1 
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Figure B.11: Angle of Attack Response of Undamaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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AoA Desired & NN Time response due to dter Plant #1 
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Figure B.12: Angle of Attack Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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Figure B.13: Pitch Rate Response of Damaged Aircraft Due to Stabilators Dipole 
Command, Single Plant, Single Network Emulation 
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qg Desired & NN Time response due to dier Plont #2 
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Figure B.14: Pitch Rate Response of Damaged Aircraft Due to Leading Edge Flaps 
Dipole Command, Single Plant, Single Network Emulation 
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q Desired & NN Time response due to dter Plont #2 
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Figure B.15: Pitch Rate Response of Damaged Aircraft Due to Trailing Edge Flaps 
Dipole Command, Single Plant, Single Network Emulation 
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nz Desired & NN Time response due to dier Plont #2 
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Figure B.16: Normal Acceleration Response of Damaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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nz Desired & NN Time response due to dter Plont #2 
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Figure B.17: Normal Acceleration Response of Damaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Single Plant, Single Network Emulation 
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AoA Desired & NN Time rasponse due to dstr Plont #2 
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Figure B.18: Angle of Attack Response of Damaged Aircraft Due to Stabilators 
Dipole Command, Single Plant, Single Network Emulation 
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AoA Desired & NN Time response due to dier Plant #2 
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Figure B.19: Angle of Attack Response of Damaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Single Plant, Single Network Emulation 
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Figure B.20: Angle of Attack Response of Damaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Single Plant, Single Network Emulation 
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Figure B.21: Spectral Response of Undamaged Aircraft Due to Leading Edge Flaps, 


Multiple-Plant, Single Network Emulation 
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Figure B.22: Spectral Response of Undamaged Aircraft Due to Trailing Edge Flaps, 
Multiple-Plant, Single Network Emulation 
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Figure B.23: Spectral Response of Damaged Aircraft Due to Leading Edge Flaps, 
Multiple-Plant, Single Network Emulation 
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Figure B.24: Spectral Response of Damaged Aircraft Due to Trailing Edge Flaps, 


Multiple-Plant, Single Network Emulation 
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Figure B.25: Pitch Rate Response of Undamaged Aircraft Due to Stabilators Dipole 
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Figure B.26: Pitch Rate Response of Undamaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.27: Pitch Rate Response of Undamaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.28: Normal Acceleration Response of Undamaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Multiple-Plant, Single Network 
Emulation 
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Figure B.29: Normal Acceleration Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Multiple-Plant, Single Network 
Emulation 
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Figure B.30: Angle of Attack Response of Undamaged Aircraft Due to Stabilators 
Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.31: Angle of Attack Response of Undamaged Aircraft Due to Leading 
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Figure B.32: Angle of Attack Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Multiple-Plant, Single Network 
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Figure B.33: Pitch Rate Response of Damaged Aircraft Due to Stabilators Dipole 
Command, Multiple-Plant, Single Network Emulation 
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Figure B.34: Pitch Rate Response of Damaged Aircraft Due to Leading Edge Flaps 
Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.35: Pitch Rate Response of Damaged Aircraft Due to Trailing Edge Flaps 
Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.36: Normal Acceleration Response Damaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.37: Normal Acceleration Response Damaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.38: Angle of Attack Response of Damaged Aircraft Due to Stabilators 
Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.39: Angle of Attack Response of Damaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 
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Figure B.40: Angle of Attack Response of Damaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Multiple-Plant, Single Network Emulation 


189 


q Desired & NN Time response 














c “LE 
plant #2: 
6 fs 8 g 10 11 eZ ae 14 Bs) 
Time [sec] Do[—-], NN[-—-] 
S 
ie 
plant #1 
6 7 8 g 10 11 12 13 14 ke 
Time [sec] Do[-], NN{--] 3 
' q Oesired NN Time response 
o 
plant Hl 
16 17 18 19 20 21 22 Zo 24 Zo 
Time [sec] Do[-], NN[--] 
q Error 
0.025 | . i 
: Ent yo Nn d\pl\\, fey. 
; pe i! We niyo oa 
| plont #2 : 3 plant #1 3 
0,02 | ! ) ! ! ! 
16 7 18 19 20 21 ZZ 25 24 ZS 
Time [sec] 0Oo[-], NN{--] 


3r 


Figure B.41: Pitch Rate Response to RBS Excitation of All Inputs, Multiple-Plant, 
Single Network Emulation 
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Figure B.42: Angle of Attack Response to RBS Excitation of All Inputs, Multiple 
Plant, Single Network Emulation 
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Figure B.43: Spectral Response of Undamaged Aircraft Due to Leading Edge Flaps, 


Multiple-Plant, Minimized Network Emulation 
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Figure B.44: Spectral Response of Undamaged Aircraft Due to Trailing Edge Flaps, 
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Figure B.46: Spectral Response of Damaged Aircraft Due to Trailing Edge Flaps, 


Multiple-Plant, Minimized Network Emulation 
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q Oesired & NN Time response due to dstr Plant #1 
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Figure B.47: Pitch Rate Response of Undamaged Aircraft Due to Stabilators Dipole 
Command, Multiple-Plant, Minimized Network Emulation 
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Figure B.48: Pitch Rate Response of Undamaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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q Oesired & NN Time response due to dter Plant #1 





fe 0:5 1 iS 2 Z.5 S) Sse) 4 4.5 =) 
Time [sec] 0O[—-], NN(-=], INunit(:] 
3 x1Q75 q Time response Errar due to dter Plont #1 


Errar 





Q oh) 1 as 2 fags. S $= 4 4.5 =, 


Time [sec] 


q Oesired & NN Time response due to dtel Plant #1 





Q 0:5 1 es) Z “5 i) 3.5 4 4.5 5 
Time [sec] 0O(-], NN{-—], iNunit{:] 
x1075 q Time response Error due to dtel Plant #1 


Error 





fe) O25 1 Lise) 2 “09 5 ro 4 4.5 5 


Time [sec] 


Figure B.49: Pitch Rate Response of Undamaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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nz Desired & NN Time response due to dler Plant #1 
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Figure B.50: Normal Acceleration Response of Undamaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 


199 


nz Desired & NN Time response due to dter Plont #1 
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Figure B.51: Normal Acceleration Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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AoA Desired & NN Time response due to dstr Plont #1 
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Figure B.52: Angle of Attack Response of Undamaged Aircraft Due to Stabilators 
Dipole Command, Multiple-Plant, Minimized Network Emulation 
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AoA Desired & NN Time response due to dier Plont #1 
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Figure B.53: Angle of Attack Response of Undamaged Aircraft Due to Leading 
Edge Flaps Dipole Command, Miultiple-Plant, Minimized Network 
Emulation 
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AoA Desired & NN Time response due to dter Plont #1 
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Figure B.54: Angle of Attack Response of Undamaged Aircraft Due to Trailing 
Edge Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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q Desired & NN Time response due to dstr Plont #2 
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Figure B.55: Pitch Rate Response of Damaged Aircraft Due to Stabilators Dipole 
Command, Multiple-Plant, Minimized Network Emulation 
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Figure B.56: Pitch Rate Response of Damaged Aircraft Due to Leading Edge Flaps 


Dipole Command, Multiple-Plant, Minimized Network Emulation 
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q Desired & NN Time response due to dter Plant #2 
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Figure B.57: Pitch Rate Response of Damaged Aircraft Due to Trailing Edge Flaps 
Dipole Command, Multiple-Plant, Minimized Network Emulation 
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nz Oesired & NN Time response due to dier Plont #2 
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Figure B.58: Normal Acceleration Response Damaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Miultiple-Plant, Minimized Network 
Emulation 


207 


oe nz ODesired & NN Time response due to dter Plant #2 


nz 





“0 0.5 1 1.5 2 2.5 S 3.5 4 4.5 


5 
Time [sec] 00[-], NN(~—], INunit([:] 
oe x1075 nz Time response Error due to dter Plant #2 
2 5 
¢ 
uw 
Q 
—5 
0 hs) 1 1.38 2 25 S 50 4 4.5 5 
Time [sec] en 
ao nz Desired & NN Time response due to dtel Plont #2 
Ps 0 
-0.2 
0 8 pes) 1 LS 2 age a) ao 4 4.5 S, 
Time [sec] 00[-], NN(-—], INunit(:] 
10 x10-5 nz Time response Error due to dtel Plant #2 
3 
Ww 





0 OS 1 ies. 2 Zo 5 hs. 4 4.5 2 


Time [sec] 


Figure B.59: Normal Acceleration Response Damaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Miultiple-Plant, Minimized Network 
Emulation 
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AoA Desired & NN Time response due to dstr Plant #2 
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Figure B.60: Angle of Attack Response of Damaged Aircraft Due to Stabilators 
Dipole Command, Multiple-Plant, Minimized Network Emulation 
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Figure B.61: Angle of Attack Response of Damaged Aircraft Due to Leading Edge 
Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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AoA Desired & NN Time response due to dter Plant #2 
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Figure B.62: Angle of Attack Response of Damaged Aircraft Due to Trailing Edge 
Flaps Dipole Command, Multiple-Plant, Minimized Network 
Emulation 
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Figure B.63: Pitch Rate Response to RBS Excitation of All Inputs, Multiple-Plant, 
Minimized Network Emulation 
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Figure B.64: Angle of Attack Response to RBS Excitation of All Inputs, Multiple- 
Plant, Minimized Network Emulation 


215 


APPENDIX C. SINGULAR VALUES PLOTS 


The graphs in Appendix C are referenced in chapter VIII. The reader ts 
referred to chapter VII for comprehensive explanation and related discussion 


regarding the plots that appear in this appendix. 


214 


Trace af singular values af weights matrix between layers 1 and 3 


Singular Values 





@) 100 200 300 400 200 600 700 800 900 eee 


# af iterations [1063] 


Figure C.1: Singular Values Plot of Weight Matrix from ComReg to comh1, Before 
Minimization 


Trace af singulor values af weights matrix between loyers 1 and 3 


Singular Values 





O 100 200 300 400 S00 600 700 800 
# af iteratians [10e3] 


Figure C.2: Singular Values Plot of Weight Matrix from ComReg to comh1, After 
Minimization 
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Figure C.3: 
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Singular Values Plot of Weight Matrix from comh1 to comh2, After 
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Troce of singular values of weights matrix between layers 2 and 5 
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Figure C.5: Singular Values Plot of Weight Matrix from OutReg to outh1, Before 
Minimization 


Trace of singular values of weights matrix between layers 2 and 5 
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Figure C.6: Singular Values Plot of Weight Matrix from OutReg to outh1, After 
Minimization 
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APPENDIX D. USERIO EMULATION PROGRAMS 


The programs in this appendix were written in C to interface with the main 
software - NeuralWorks Professional H/PLUS™ by NeuralWare, Inc. that used in this 
research. There are seven programs involved here. The main program emul.c is the 
USERIO program and is based on the general network structure described in Chapter 
VI. There are five dedicated header-files included. 

The first is matlab_link.h which makes the interface with MATLAB™ as a 
server program to calculate the singular values of the weight matrices while the main 
program is running. The second header file is task depended. The task shown here 
is emulation. The F/A-18A undamaged and damaged plant that were emulated are 
simulated in a user defined header-file /78/_modelN.h. This file may be edited for 
different plants with the restriction described in the program listing. 

The third header-file is net_paramN.h includes all the user accessed parameters 
that influence the network structure and task dependent calculation. The forth header- 
file net_util.h defines parameters and functions necessary for the main program. The 
last file, window_dat.h defines parameter and user graphic routines. All the 
parameters and variables are explained in the body of the program mostly with self- 
explanatory names. 

The control strategy program emulbpt.nnc controls the flow of the main 
USERIO program, and applies the back-propagation algorithm. The use of emul.c 


must be in conjunction with this control strategy. 
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re ee eg ee ee ee ee A eR RRA RRR RRR 


Ae $Author : $ Shahar Dror =] 
= $Date : $ 22 Sep 1991 */ 
/* $Source : $ emul.c sald 
ifs $Revision : $ 16 Jun 1991 */ 
l= The NN is trained to produce the next plant state based on */ 
ih the current and past input commands and real-plant outputs */ 
ve Whenever adding/changing Tasks, search for ‘Task’ to a7 

ibe accommodate/verify correct interpretation! Search also all */ 
7* #include files as well!! */ 


aK a ha Ree eS eA ee eee ee A eee eX / 


#include <stdio.h> 
#include <math.h> 
#include <sys/types.h> 
#include <string.h> 
#include <malloc.h> 


/* Definitions */ 


/* sz=ssss=s===2 */] 

#def ine UIO_SERVER 1 
#define SERVER_EMULATOR_FILE 1 
#define SKIP COMPAT DATA 1 


#define BEEP (printf(''*G"), fflush(stdout) ) /* Sound a Beep */ 


#define SET 1 

#define RESET 0 

#define LEARN 0 /* LEARN mode indicator */ 

#define TEST 1 /* TEST mode indicator */ 

#define IClearn 2 /* Initial conditions in LEARN mode indicator */ 
#define ICtest 3 /* Initial conditions in TEST mode indicator */ 
#define TEST1 4 /* TEST mode, input type is RBS indicator */ 
#define JACOB 5 /* Jacobian calculation identifier */ 

#define In_lay 0 /* Input layer */ 

#define ComReg_lay 1 /* Command(input) regression layer */ 

#define OutReg lay 2 /* Network output regression layer */ 

#define classification 5 /* Define task mode */ 

#define emulation 1 /* Define task mode */ 

#define combined_cl_em 2 /* Define task mode */ 

#define controlling 3 /* Define task mode */ 

#define cntriThruPlant 4 /* Define task mode */ 

/*#define Task classification /* Define operating task mode */ 
#def ine Task emulation /* Define operating task mode */ 


#include "“userutl.h" 
#include "uio_stre.h" 
#include "uio_tran.h" 
#include “ulio reqs.pro" 
#include “uio_cmds.pro"” 
#include "matlab link.h" 


#if Task==emulation /* Set up emulation mode */ 
#include "f18l_modetN.h" 

#endi f 

#if Task==controlling /* Set up controlling mode */ 
#include "ref_modetN.h" 

#endi f 


#include "net_paramN.h" 
#include "net_util.h" 
#include "window _dat.h" 


#define MAXRAND (Ox4fffffff1) 
#define Out_lay (nlayp-1) /* Output layer number */ 
#define OutBuf_lay (Out_lay-1)  /* Output buffer layer number */ 


/* Declarations */ 


/* sssssssesesse= */ 

extern long random(); /* Random number generator */ 

int Leycle=0,Rcycle=0; /* Learning, Recall/Testing cycles counter, */ 

int LR_flag; /* Learn(1)/Recall(0) flag */ 

char *sp; /* String pointers */ 

time_t *tloc; /* Pointer to current time for random # initialization */ 
ihe 
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[REAERRRERARRAREREREREREAERAEREE EAD RA AO ee 


/* Functions necessary for handling the User 1/0 package itself. */ 
[RRRRERERARAE RARER RES HERS SESS Se ee 


NINT UIO Init(file) 
TEXT *file; 


NINT ret_val = UIO_OK; 
int mnr; /* Model number index */ 
long ctrp(8); /* An array of 8 to interact with control-strategy*/ 


srandom(time(tloc)); /* Random # initialization, using a new seed each start */ 
/* Initialize counters which interact with control strategy */ 
ug_ctrio(ctrp,0);ctrp(7]=ctrp(6)]=ctrp(4]=0;ug_ctrio(ctrp, 1); /*aux2d=wrk 1=wrk2=0*/ 


if (check_net()==UIO_ERROR) return(UIO_ERROR); /* Network validation */ 

if(start_matlab()==UIO_ERROR) return(UIO_ERROR); /* Linking to matlab */ 

for(mnr=1;mnr<=Nmodel s;mnr++) /*Loading models matrices*/ 
if (load_model (mnr)==UIO_ERROR) return(UIO_ERROR); 


/* Load system model if exists */ 
ifCeffector_lay_length>0)¢ 
for(mnr=T;mnr<=Nref_models;mnr++){ 
if (load_ref_model (mnr)==UIO_ERROR) return(UIO_ERROR); 
} 


set_net(); /* Set and initialize network and model related parameters */ 
ret_val=UIO_Attention(); /* Establishing working format */ 
if(ret_val==UIO_ERROR) PutStr("Exiting... Restart the run."); 


return(ret_val); 
} 


NINT UIO_Term(process) 
NINT process; 


{ 

close_files(); /* Close all opened files */ 

stop_matlab(); /* Unlinking matlab and closing relevant files and pipes */ 
if(LR_flag==LEARN) PutFmt("\nTraining cycles: %i",Lcycle); 
if(LR_ flag==TEST ) PutFmt("\nTest cycles: %i",Rcycle); 

PutStr( " bye bye\n" ); 

) 


NINT UIO Attention() 


{ 

NINT ret_val = UIO_OK; 
int tmin, tmax; /* Min & max time levels for H->L swept square wave */ 
int last_LR; /* The value of LR_flag when entering UIO Attention. */ 


intrrpt_flag=SET; /* Setting interrupt flag. Causes network checking */ 
close_files(); /* If files were open, close them */ 
/* Inquire mode type LEARN/TEST */ 
last_LR=LR_flag; 
LR_again: 
PutFmt("\nIs this run a LEARN(O) or TEST(1) <#!=0,1 - quit> ? (%i] ",LR_flag); 
sp=GetStr(); 
if(*sp) sscanf(sp,"%i",&LR_ flag); 
if(!C(LR_flag==0 || LR_flag==1))¢ 
PutStr("\nLR_flag must be 0 or 1. Do you want to quit? [y]/n "); 
sp=GetStr(); 
LR_flag=last_LR; 
if(!*sp ji, *sp==’y’) return(UIO_ERROR); 
else goto LR_again; 


) 


model _nbr=1; /* start training with unimpaired (first) model */ 
rand_models (0] =0; 
if (Nmodels>1){ 
model_nbr=0; 
Againe: 
PutFmt("\nEnter Model #4 to train/test (1-%1), O-for all (%i] ", 
Nmodels,model_nbr); 
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sp=GetStr(); 
if(*sp) sscanf( sp, "%d"", &model_nbr); 
if( model_nbr >Nmodels || model_nbr<0O )¢ 
PutFmt(" %s",input_name[0] ); 
goto Again2; 


} 
Almodel_flag=((model_nbr!=0)?SET:RESET); 
de , 
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if (LR_flag==TEST)¢ /* TES) mode */ 


/* User select input to be used */ 
Input_Type_again: 
PutFmt("\nEnter Input Type 1) %s, 2) %s, 3) %s", 
input_name [1], input_name [2] , input_name [3] ); 
PutFmt("\n4) %s, 5) 4s, 6) %s ([4i] ", 
input_name [4] , input_name [5] , input_name [6] , input); 
sp=GetStr(); 

if(*sp) sscanf( sp, "%d", &input); 

if( input >Ninput_types ;{ input<1 )¢ 

PutFmt(" %s",input_name(Q) ); 

goto Input_Type_again; 


) 
PutFmt("\nInput: %s selected", input_name([input] );/*Display selection*/ 


netSave flag=RESET; /* Resetting best network saving flag */ 


Rcycle=0; /* Initializing test cycles counter */ 
min_RMS=1.0; max_RMS=0.0; /* Initializing rel_err extrimums */ 
if ( input==1)¢€ /* Random Binary Sequence */ 


if( Nmodels>1 && model_nbr==0)¢ 
if(Nmodels>2) rand_array(Nmodels, rand_models); 


else{ 
rand_models[(Q0)=2; rand_models(1]=1; 
) 
PutFmt("\nEnter max random cycles to switch models (%i]) ",rmRC); 
sp=GetStr(); 
if(*sp) sscanf(sp,"%i", &rmRC); 
mRC=rmRC; 


} 
else mRC=Large number;/* Select a large # */ 
inTest=0; /* Test on all inputs on all inputs */ 


) 
if(model_nbr==0) model_nbr=1;/* Start with original model (for test) */ 


ifCinput==2 |} input==3 |} input==4 || input==5)¢ 
/*Sequential Random Binary, Swept square/sine wave inputs */ 
PutFmt("\nEnter upper TEST frequency (%f Hz] ", Fmax); sp=GetStr(); 
if(*sp) sscanf(sp,"4lf",&Fmax); 
if(Fmax> 0.5/ts) 
PutFmt("\nUpper freq exceeded. Fmax set to %f Hz'', Fmax=0.5/ts); 
PutFmt("\nEnter lowest TEST frequency (%f Hz) ",Fmin); sp=GetStr(); 
if(*sp) sscanf(sp,"4lf",&Fmin); 
mRC=( int )Cexpe(ceil(-logé(ts*Fmin)))); 
PutFmt("\nSet Test cycles >= 4li", 
mRC*forced_excitation_length*(Nmodels! =0?Nmodels:1)); 
if (input==3)¢ 
tmin=floor(0.5/Fmax/ts); 
tmax=ceil(0.5/Fmin/ts); 
if (tmax>mRC*0.15) tmax=ceil(mRC*0.15); 
if(sweep_HL(tmin, tmax,mRC,sweep_change_times )==UIO_ ERROR) 
return(UIO_ERROR); 
Sweep_count=1;Nsweep_change=0; /* initializing Sweep counters*/ 
) 
ifCinput==4 |} input==5)¢ 
dmin=log10(Fmin); 
df log=( log10( Fmax)-dmin)/(mRC-1); 
df | in=( Fmax-Fmin)/(mRC- 1); 
) 
inTest=0; 
) 


if (input ==6) ¢ /* One second dipole input */ 
mRC=(int)(dpl (3) /ts); 
PutFmt("\nSet Test cycles >= 4Li", 

mRC* forced_excitation_length*(Nmodels!=0?Nmodels:1)); 


inTest=0; 
) 

PutStr("\nDo you want to save TEST results? [y]/n "); sp=GetStr(); 
if(!*sp |! *sp=='y/ )¢ /* open files for test results */ 


if(open_TEST_files()==UIO_ERROR) return(UIO_ERROR); 
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else 


— a oe ee ee 
—_ oe oe oe oe oe oe oe 


if(InPut && NNo && DESo) TestDataSave_flag=SET; 


TestDataSave_flag=RESET; 


/* Don’t save TEST results */ 
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=s====\/ ==Start====ss==LEARNssassese5\ /sesesss=_LEARN=====start=====\/===* 
Leycle=0; 
input=1; /* Training with Random Binary Sequence */ 


PutFmt("\nInput: %s selected", input_name(input] );/*Display selection*/ 
if( Nmodels>1 && model _nbr==0)¢ 
if (Nmodels>2) rand_array(Nmodels,rand_models); 
else{ 
rand_models(0]=2; rand_models(1]=1; 
} 
Switch_again: 
PutFmt("\nEnter min random cycles to switch models (%1] ",minSwW); 


sp=GetStr(); 
if(*sp) sscanf(sp,"%i", &minSW); 
maxSW-=1; 
PutFmt("\nEnter max random cycles to switch models (%i] ",maxSw); 
sp=GetStr(); 
if(*sp) sscanf(sp,"%i", &maxSw); 
maxSW+=1; /* In order to reach maxSW value */ 


i f ( (maxSW-minSW)<=0)¢ 
PutFmt(" %s",input_name(0] ); 
goto Switch_again; 
2 
mRC=mi nSW+random( )4(maxSW-minSW); 
model_nbr=1; 
> 
else{ 
mR C=maxSW=Large_number-2;minSW=maxSw- 1; 
} 
PutStr("\nDo you want to use Reset-Training mode ? [y]/n "); 
sp=GetStr(); 
if(!*sp {| *sp==’y’)¢ 
ResetTrain_flag=SET ; 
PutFmt("\nEnter Reset-rate (in learn-cycles) per model (%1]: ", 
ResetTC); sp=GetStr(); 
if(*sp) sscanf(sp,"%i", &ResetTC); 
nRTC=ResetTC*(Almodel_flag==SET?1:Nmodels); 
) 
else ResetTrain_flag=RESET; 
PutStr("\nDo you want to Reset if output exceed limits ? [y]/n "); 
sp=GetStr(); 
if(!*sp jj, *sp==’y’) ResetLimit_flag=SET; 
else ResetLimit_flag=RESET; 


if€forced_excitation_length>1){ 
Input_number_again: 
PutFmt("\nEnter Input # to train (1-41), O-for all (%i) " 


,forced_excitation_length, inTest); sp=GetStr(); 
if(*sp) sscanf( sp, "%d", &inTest); 
if¢ inTest >forced_excitation_length || inTest<0O )¢ 


PutFmt(" %s",input_name(0] ); 
goto [nput_number_again; 
2 


) 
PutFmt("\nEnter file name to save SVD (%s] ",svdn); 
sp=GetStr(); 
if(*sp) sscanf(sp,"%s",svdn); 
fputs("clear \n",matlab); /* Clearing matlab working space */ 
netSave_flag=RESET; /* Initially don’t save best network */ 
PutStr("\nDo you want to save best network ? (yJ/n "); sp=GetStr(); 
if(!*sp |] *sp=='y!)¢ 
netSave_flag=SET; 
PutFmt("\nEnter convergence threshold (4lf]: ", convergence); 
sp=GetStr(); 
if(*sp) sscanf(sp,"4lf", &convergence); 
) 


min_sum=convergence*NErrConv; /* Initialize best-net saving */ 
Rstart_flag=RESET; /* Didn’t go through UIO Recall start (yet). */ 


/*ssssssses/ \====end=========_LEARN=sssessas=/\=ssses==_L EARN==s==send===s===/\===*/ 
PutStr("\n UIO_Attention finished. You may continue... "); 


} 


return(ret_val); 
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Jf TAEAKRARREAKKKKAKAKKARNKNR WAKA RKKKHRKKAKKK CHAKA KERR RKE RAR REK HEH / 


/* Functions necessary for handling a learning session. hall f 
J NARA AAAAAARAAKKKKKRRRK RARE WEE belied aeladialiahalia halal Ralele Rahal Ra lelelshahahehehahehehehehelalehelel | 


NINT UIO Learn _Start() 
/* Learn start; occurs once at each start (also after intermediate interruption. 
Prepares data for next coming cycle */ 


{ 
NINT ret_val = UIO_OK; 
int hist,hist_in, hist _out; /* input_data vector initialization counters */ 
int mnr; /* Model number counter */ 
int 1,),K3 /* General indices */ 
long ctrp(8]; /* An array of 8 to interact with control-strategy*/ 


error_flag=RESET; /* Resetting error flag. No error yet! */ 


/* Check compatibility of initialization and session type (L/R) */ 
if(LR_flag!=LEARN && Rstart_flag==RESET){ /* Wrong. exit! */ 
PutStr("\n!! You run a LEARN session, but initialized for RECALL."); 
BEEP; /* sound a warning Beep */ 

PutStr("\nFor LEARN, run I/0...>Attention I/O. For RECALL, run Run. .>Recall"); 
PutStr¢(". Press <RETURN>.'); 
GetStr(); return(UIO_EXIT); 
> 
ifCintrrpt_flag==SET)¢ /* Validate network compatibility */ 
if(check_net()==UIO_ ERROR) return(UIO_ERROR); 
intrrpt_flag=RESET; 
> 


if(LR_flag==LEARN)¢ /* Learning mode */ 
if(Lceycle%mRC==0 || netReset_flag==SET){¢ 
for(mnr=1;mnr<=Nmodels;mnr++){ po boad Cn */ 


desired | output (IClearn, mnr,excite,&desout (mnr-1) [0] ); 
/* Check for errors in desired _output¢) oi 
ifCerror_flag==SET) return(UIO_ERROR); 
ifceffector _lay_length>0){/* Initial Effector ball / 
for(i=0;i<effector_lay_length; i++) 
Effector tmr- 1] [iJ =0.0; 


y 


if(LR_flag==TEST){ /* Testing mode */ 
Tf(Reycles“mRC==0 && input!=1)¢ 
ifCinTest>=forced_excitation_length){ 
if((Almodel_flag==SET |! ++model_nbr>Nmodels) && 
UIO _Attention( )==UI0_ERROR) return(UIO_ ERROR); 
else inTest=0; 


> 
/*Load test IC*/ 
desired_output(ICtest,model_nbr,excite, &desout [model_nbr-1]} (0]); 
/* Check for errors in desired | output () Ne) 
ifCerror_flag==SET) return(UIO ERROR); 
if(effector_ lay_length>0)¢ i> Initialize Ef fector*/ 
for(1=0;i<effector _lay_length; i++) 
Effector [model _nbr-1] [7] =0.0; 


> 

ifCinput==3) ¢ /* initializing sweep counters*/ 
sweep_count=1; Nsweep_change=0; 

> 


ifCinput==4 || input==5)¢ 
ff log=dmin; /* Initializing (log10) frequency */ 


flin=Fmin; /* Initializing linear frequency */ 
T=0; /* Initializing Time */ 
ifCinput==6) T=0; /* Initializing Time */ 


y 

/* Place the initial (measurements) conditions in the input_data vector */ 
if((LR_flag==TEST && Rcycle%mRC==0) |! (LR_flag==LEARN && 
(Leycle“mRC==0 |! netReset_flag==SET)))¢ 


for(mnr=((LR_flag==LEARN && Almodel_flag==RESET)?1:model_nbr); 
mnr<=((LR_flag==LEARN&&A1model_flag==RESET)?Nmodels: model_nbr); mnr++){ 


Us p255) 


/* Set counters for cases in which input_history!=output_history */ 
if(in_history>out_history){ 

hist=in _history; “hist_ in=0; hist_out=out_history-in_history; 

) 

else¢ 

hist=out_history; hist_out=0; hist_in=in_history-out_history; 


for( j=0; j<=hist;j++)€ /* Fill in input and regression layers*/ 
/* Last pass @ j=hist is preparing for first presentation, not loaded to REG 
UCO) and Y(O) Cinitial values/conditions) are lost */ 


/* Create excitations */ 

for(++hist_in, 1=0; i<forced_excitation_length; i++)¢ 
if (Tnput= =1 && (j%RbsUd)==0)/* Random Binary Sequence */ 
excite(i]=(inTest==0?RBS:(inTest==(1+1)?7RBS:0.0)); 
ifCinput==2 && ( j~ARbsUd)==0) 
excite({i]=(inTest==i?RBS:0.0); /*Seq RBS*/ 
if (input==3) excite(i]=(CinTest==i?(one=1.0):0.0);/*HL*/ 
ifCinput==4) excite(i]=CinTest==171.0:0.0);/*SwL-HSqr*/ 
if(input==5) excite([i]=0.0; /*Swept L-H Sine Wave */ 
if(input==6) excite{i]=0.0; /*Dipole */ 
/* To induce real zero, subtract the bias. */ 
ifCexcite(i]==0.0) exciteli]=-bias_I(mnr-1) (il; 

input_data(mnr-1] (i]=exciteli]; 


> 
if(unforced_excitation_length>0)¢ 
unforced |_excitation(LR _ftag,mnr,&input_data(mnr-1] (0), 
& input_data (mnr- 1] [forced_ excitation length] ); 
/* Check for errors in unforced_excitation() ) */ 
if(error_flag==SET) return(UIO_ERROR); 
> 


if(type_length>0){ /* Set to Active the model in tact */ 
Type((LR_flag==LEARN?1Clearn: [Ctest),mnr,&desout (mnr-1] (0), 
&E ffector (mnr-1] (0) ,&input_data(mnr-1] [excitation_length] ); 
/* Check for errors in Type() ) */ 
if(error_flag==SET) return(UIO_ERROR); 


/* Now load into ComReg the relevant inputs */ 
ifChist_in>O && j<hist)¢ 
for(i=0; i<regressed_ input_ length; i1++)¢€ 
k=(hist_in- 1)*regressed input_length+i; 
ComReg (mnr- 1] (k]=input_data(mnr- 1] [RIsort lil]; 
> 


> 
/* Put previos output values int OutReg */ 
if(++hist_out>0)¢ 
for(i=0;i<regressed _output_length; i++)¢ 

if(j<hist)€ 
K=1+(hist_out-1)*regressed_output_length; 
OutReg [mnr -1] [k] =desout (mnr-1] (ROsortli]]; 
> 


else 
desoutOld(mnr-1) (1) =desout (mnr-1) (ROsort(i]]; 
) 


/* Using current excitation, calculate desired_output at t=(current+1)*/ 
desired_output(LR_flag,mnr ,&input_data(mnr-1] (0] ,&desout (mnr- 1] (0) ); 
> /* (end for-j) Do until ‘both regression vectors are filled */ 

> /* (end for-mnr) Do for all relevant models */ 
ug _ctrio(ctrp,0); /* read counters content */ 


ctrp(6]=1; /* Set wrki=1 to load Regressors values */ 
ug_ctrio(ctrp,1); /* write counters content */ 
} /* (end if-%mRC==0) */ 

netReset_flag=RESET; /* Do not reset network next time */ 


/* Display the starting condition */ 
1f(LR_flag==LEARN) PutFmt("\nTraining cycles: %4i",Lcycle); 
if(LR_ _flag==TEST) PutFmt("\nTesting cycles: %i",Rcycle); 
PutFmt(", Input type: %s, model #: %i, input # in effect: %i", 
input_nameCinput], ‘model_nbr, (input==1?inTest:(inTest+1))); 
return(ret_val); 
) 
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NINT UIO_Learn_Input(LayN, nPEs, Datap) 
NINT LayN; 

NINT nPEs; 

SREAL *Datap; 


NINT ret_val = UIO OK; 
int i; 7* Counter index */ 
/* WARNING: Layers that are not to be filled, are to be controlled via the 
control-strategy, and not by not filling them in UIO_Learn _Input, since 
this function will write the last data that was addressed in this function */ 


/* Fill the input layer with the required input data vector */ 
if(LayN==In tay) 
for(1=0;i<nPEs;Datap(iJ=input_data[model_nbr-1] [iJ], i++); 


} 
/* Fill the Command Regression layer with the last output(model#) data vector */ 

if(LayN==ComReg_lay && regressed_input_length>0){¢ 
if(regressed_input_history_length<nPEs){ 
PutFmt("\n ComReg vector too small. Vector length :%i, nPEs=%i"', 
regressed_input_history_length,nPEs); 
PutStr("\nPress RETURN to stop and exit"); GetStr(); 
return(UIO_ERROR); 
> 


for(i=0;i1<nPEs;Datap(i] =ComReg[model_nbr-1)(i], i++); 


/* Fill the Output Regression layer with the last output(model#) data vector */ 
if(LayN==OutReg lay && regressed_output_length>0){ 
if(regressed | output_ history_ length<nPEs){ 
PutFmt("\n OutReg vector too small. Vector length 241, NPEs=Ai", 
regressed_output_history_length,nPEs); 
PutStr("\nPress RETURN to stop and exit"); GetStr(); 
return(UIO_ERROR); 
> 
for(i=0;1<nPEs;Datap[i]=OutReg([model_nbr-1] Ci], i++); 


} 
/* Fill the Output Buffer layer with last desired output(model#) data vector */ 
if(LayN==OutBuf lay && regressed_output_length>0){ 
for(i=0; 1<nPEs;Datap[i]=desoutOld[model_nbr-1] Ci], i++); 


} 
return(ret_val); 


NINT UIO Learn_Output(LayN, nPEs, Datap) 
NINT LayN; 

NINT NMPEs; 

SREAL *“Datap; 

{ 


NINT ret_val = UIO OK; 
int i; /* Eouncer index */ 


/* Fill output layer with the desired output data (sensors measurements) */ 
for(1=0; 1<nPEs; i++) 
Datapli]=desoutOld{model_ nbr- 1) [i] =desout [model_nbr-1] [i]; 
return(ret_val); 


/* 


yoped| 


NINT UIO Learn_Result(LayN, mPEs, Datap) 


NINT LayN; 

NINT mPEs; 

SREAL *Datap; ibe: Datap = pointer to array that contains output data */ 
NINT ret_val = UIO_OK; 

double rel_err,err_MS; /* Relative output error */ 

double abs_msr; /* Absolute value of the measurement */ 

Static int ~cycl (Nmodels+1], iInput_index=0,model_index=0; /* Cycles counter */ 

static int his _cycl [2] (Nmodels+17; Pes History cycle counter for Learn and Test eld 

long ctrp[8]; /* A pointer to array of 8 counters that control NW run */ 

int last “model: /* In TEST, & input =1, saves previos model number */ 

Int 1,)- /* Counter index */ 


double nn[output_lay_length] ; 


/* This function reads the output of the current step, 
and prepares the input for the next step */ 


ifCerror_flag==SET) return(UIO_ERROR); 

if(effector_ lay_length>0 && LayN==effector_lay){/*Get control outputs*/ 
if(mPEs!=effector_lay_length)<{ /* Check matching */ 
PutStr("\nWrong control-layer or control-length mismatch"); 
PutFmt(". LayN=%i1, nPEs=41, effector-length=%i", 

LayN,nPEs, effector _'ay_ Length); 
efLay_flag==SET; /* Enable setting Effector _lay” # */ 
check_net(); /* Resetting effector_lay # */ 
PutStr("Rerun after termination. "); 
return(UIO_ EXIT); 

/*end if(nPEs...*/ 
for(i=0;1<nPEs;Effector [model_nbr-1] Ci] =Datap[li], i++); 
/* if effector layer is the output layer, continue */ 
if¢effector_lay != Out_lay) return(ret_val); 
} /*end ifceffector...*/~ 
/* Increment the counter and update the display as necessary */ 
if(LR_flag==LEARN){ /* Display Learn cycles in LEARN mode */ 
if(ResetLimit_flag==SET && LayN==Out_lay && 
Limits(model_nbr ,Datap)==0)¢ 


ug_ctrio(ctrp,0); /* read counters content */ 
ctrp(7)=1; /* Set wrk2=1 to reset network */ 
ug _etriocetrp, 1); /* write counters content */ 


PutStr("\nOutput Limit Exceeded. Resetting..."); 
netReset_flag=SET; 
return(UIO_ Learn | Start()): 


. 
i f(LayN==ComReg_lay){ /* Save ComReg(model#) output for next*/ 
if(ResetTrain_flag==SET && Leycle%nRTC==0)¢ 
netReset_flag=SET; 
return(UIO_Learn_Start()); /* Reset plant/ref */ 
) 


for(i=0; 1<nPEs;ComReg (model_nbr-1) (i) =Datap(i] , i++); 
return(ret_val); 
) : 
/* Save OutReg(model#) output for next*/ 
if(LayN==OutReg lay){ 

if(ResetTrain_flag==SET && Lcycle%nRTC==0) 

return(ret_val); 
if(regressed output_length>0) 
for(i=0;1<nPEs;OutReg (model _nor-1) (i) =Datap(i), i++); 

return(ret_val); 


Lcycle++; 
if ((Leycle%500)<1) PutFmt("\nTraining cycles: %i",Lceycle); 


) 
if(LR_flag==TEST)¢ ° 
if(TestDataSave_flag==SET){ /* Save data for analyzing results*/ 
/*Save input commands in TEST mode*/ 
ifCinput==1) /* Store all inputs */ 
for(i=0;i<forced_excitation_length; i++) 
fprintf(InPut, "% 1.5e ", input_ data (model _nbr-1) (1) ); 
else 
fprintf¢(InPut,"% 1.5e ", input_data(model_nbr-1] (inTest] ); 
fprintfCInPut,"\n"); 
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/*Save desired-outputs in TEST mode*/ 
for(i1=0; 1<noutp; i++) 
fprintf(DESo,"% 1.5e ",desout [model _nbr-1] [i] ); 
fprintf(DESo,"\n"); 
/* Save the NN outputs */ 
for(i=0;1<nPEs;fprintf(NNo,"% 1.5e ",Datap[i] ), i++); 
fprintf(NNo,"\n"'); 


) 
ug_ctrio(ctrp,0); /* read couters content */ 
ctrp[6] =0; /* Set wrk1=0. Do not load Regressors values */ 
ug _ctrio(ctrp,1); °/* write couters content */ 
Rcyclet+; 


if(Rceycle%”mRC==0) ¢ 
if(input==1 && Nmodels>1)¢ 
last_model=model_nbr; 
model_nbr=rand _ models (model index++]; 
PutFmt("\nTesting cycles: i" »Reycle); 
PutFmt(", model #: %i1, input # in effect: BAe 
model_nbr, (input==1?inTest:(inTest+1))); 
/* Prepare values for OutBuf layer 
(for test with desired values!) */ 
for(1=0;i<output_lay_length; i++) 
desoutOld{model _nbr- 1) Li] =desout [last_ model-1] Cil; 
if((model_index=model_index%Nmodels)==0 && Nmodels>2) 
rand_array(Nmodels,rand_models);/*new rand modelseq*/ 


else{ 

inTest++; 
return(UIO_Learn_Start()); 
) 


> 
if((Rcycle%500)<1) PutFmt("\nTest cycles: %i",Rcycle); 
2 
cycl (model_nbr]++; /* Increment excitation update cycle counter */ 
/* Generate new control- input and place it in the first elements of input_data*/ 
if ( input==1){€ /* Random binary sequence */ 
if(Ccycl [model _nbr]=cycl (model _nbr] ZRbsUd)==0)¢ 
for(i=0; -i<forced_ excitation_length; i++){¢ 
excite[i]= (inTest==0?RBS: (inTest= =|+12RBS: 0.0)); 
/* To induce real zero, subtract the bias. */ 
if (excite([i] ==0.0) excitelil=-bias_I [model_nbr-1] (i); 
input_data[model_nbr-1] (i]=excite[i]; 
> 
> 


) 

ifCinput==2 && (cycl [model_nbr]=cycl [model_nbr] %RbsUd)==0) 
excite(linTest] =RBS; 

if( input==3) exciteLinTest] =SweptSquareWaveH2L(); 

ifCinput==4) excite[inTest] =SweptSquareWavel2H(Rcycle); 

ifCinput==5) excitelinTest] =SweptS ineWaveL2H(Rcycle); 

if (input==6)¢ 


T+=ts; 
excitelinTest]=DIPOLE((T)); /* 1 second dipole */ 
> 

ifCinput!=1)¢ /* Update input data with the new control */ 


/* To induce real Zero, subtract the bias. */ 
if(excitelLinTest] = =0. 0) excitel[inTest]=-bias_I [model_nbr-1] [inTest]; 

input_data(model_nbr-1] [inTest]=excitelinTest]; 
> 
/* If input consist of unforced excitation, calculate it. td) 
if (unforced_excitation_length>0) 
unforced_excitation(LR_flag,model_nbr,&input_data(model_nbr-1] [0], 

&input_ data(model_nbr- 1] [forced_ excitation length] ); 
/* If input consist of type indication, calculate it. “*/ 
if(type_length>0){ /* Set to Active the model In tact */ 

or(i=0;1<nPEs;nn[i]=Datap[i] , i++);/*convert float to double */ 
Type(LR_flag, model _nbr, nn 8 ffector [model —nbr-1] (0), 
& input _ dat a({model _nbr- 1] [exci tation. length] ); 


/* Save best network if desired (convergence test) */ 
if(((netSave_flag==SET |} (Nmodels>1 && rand_models [0] !=0))&&inTest==0) 
&& LR_flag= =LEARN && Lcycle>NErrConv){ 
for(i=0,err_MS=0.0; i<noutp; i++){ 


Lo 


rel_err=(desout [model_nbr- 1] [i] -(double)Datap[i] ); 
err_MSt+=rel_err*rel_err; 


) 
if(Cerr_RMS=sart(err_MS))<=convergence){ 
i f(++conv _test [0] ==Lcycle)¢ 
for(i=NErrConv; i>1;sum[i]=sum[i-1],i--); 
sum(0]+=(sum{1]=err_ RMS); 
1f(++conv_test [1] >=NErrConv){ 
if¢sum([0]<min_sum && netSave_flag==SET)¢ 
SaveNet (BestNetName, NetSaveMode) ; 
min_sum=sum(0] ; 
PutFmt("\nTraining cycles: 41",Lcycle); 
PutFmt(", Input type: %s, model #: Zi", 
input_name [input] ,model_nbr); 
printf("\nTraining cycles: %4i",Lcycle); 
printf(", model #: %1, mean_min_RMS=%f", 
model _nbr,min_sum/NErrConv); 
fflush(stdout); 


> 
sum[(Q) -=sum(NErrConv] ; 


else€ 

conv_test (0) =Lcycle; 
conv test[1]= 1 
sum(0]=err _RMS; 

) 


) 


/* Prepare output for next presentation of that model */ 
desired_output(( input==1&&LR_flag==TEST )?TEST1:LEARN,model_nbr, 
&input_data [model _nbr-1) (0), &desout [model_nbr- 1) (0) ); 
/* Switch models Tf it is the time to do so. */ 
ifCinput==1 && LR flag==LEARN && Lceycle%mRC==0 && Nmodels>1)¢ 
model_nbr=rand_models [model _ index] ; 
if ((model _index=((++mode | _index)%Nmodels))= =0 && Nmodels>2) 
rand |_array(Nmodels, rand | models); /*new random model -seq*/ 
mRC=L3cyc l e+minSW+random( )%(maxSw- minSW); 


return(ret_val); 
) 


NINT UIO_Learn_End() 
{ 
NINT ret_val = UIO_OK; 
if(LR_flag==LEARN){ 
/* Save matlab workspace containing SVOs matrices */ 
fputs("save ",matlab);fputs( svdn ,matlab); fputs("\n",matlab);/* save */ 
PutFmt("\nTraining cycles: %i",Lcycle); 
> 
if(LR_flag==TEST) PutFmt("\nTest cycles: “%i",Rcycle); 
intrrpt_flag=SET; /* An interrupt occurred. Check network. */ 
PutStr(". To terminate press <RETURN> "); 
GetStr();PutStr¢("\b ... Done! ")>; 
return(ret_val); 
> 
7* 
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/* Functions necessary for handling a recall or testing session. */ 
[RARRRREREREREREREKKKEKK HEKKKKKRKREKREKRERERER ERK KR EKERERERRERERERERER f 


NINT UIO Recall Start() 
€ 


NINT ret_val = UIO_UPDATE; 
/* Check compatibility of initialization and session type (L/R) */ 


if(LR_flag!=TEST)¢ /* Wrong. Exit! */ 
PutStr("\n!! You run a RECALL session, but initialized for LEARN."); 
BEEP; /* sound a warning Beep */ 


PutStr("\nFor RECALL, run I/0...>Attention I/0. For LEARN, run Run...>Learn"); 
PutStr(". Press <RETURN>."); 
GetStr(); return(UIO EXIT); 
> 
Rstart_flag=SET; 
return(UIO Learn_Start()); 
} 


NINT UIO Read(LayN, nPEs, Datap) 
NINT LayN; 

NINT nPEs; 

SREAL *Datap; 


{ 
/* In TEST mode load desired outputs */ 
if(LayN==Out_lay) return(UIO_Learn_Output(LayN, nmPEs, Datap)); 
/* In TEST/Recall mode get last values */ 

‘ else return( UIO_Learn_Input(LayN, mPEs, Datap)); 


NINT UIO_Write(LayN, nPEs, Datap) 
NINT LayN; 

NINT nPEs; 

SREAL *Datap; 

{ 


NINT ret_val = UIO_OK; 
return( UIO_ Learn_Result(LayN, nPEs, Datap)); 


NINT UIO Write Step(LayN, nPEs, Datap) /* Not used in BackPropagation */ 
NINT LayN; 

NINT mPEs; 

SREAL *Datap; 

{ 


NINT ret_val = UIO OK; 
return(ret_val); 


NINT UIO_ Recall _Test(LayN, nPEs, Datap) 
NINT LayN; 

NINT mPEs; 

SREAL *Datap; 


{ 
return(UIO Learn_Output(LayN, mPEs, Datap)); 
} 
NINT UIO Recall _End¢) 
€ 


NINT ret_val = UIO_OK; 
Rstart_flag=RESET; 

return(UIO_Learn_End()); 
, 
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[ARRAS RAAAAR EERE ne ee ee eed 


/* Other miscellaneous functions. ay 
[ERAERERERERERAEARESRAEAN AN ERR EN AA Ee ee ee 


NINT UIO_Instrument(Instrument_id, nOataElems, DataElemp) 

NINT Instrument_id; 

NINT nDdataElems; 

SREAL *DataElemp; 

{ 

NINT ret_val = UIO OK; 

double dDE(nInstElmnt]; /* an array to convert DataElemp from float to double */ 
char inst(maxInstName],c; /* Variable of an instrument name */ 


int 1; /* Counter index */ 

int rows,cols; /* # rows, # columns in weight matrix */ 
if(mDataElems>ninstElmnt){ /* Check d@DE is big enough */ 
PutFmt( 


"\nInstrument transfers too many points. ninstElmnt=%1, but has to be:%i" 
ninstElmnt,nDataElems); 

printf 

"\nIinstrument transfers too many points. niInstElmnt=4i1,but has to be:%i" 
,ninstelmnt, nataElems); 

error_flag=SET; PutStr("\nPress RETURN to stop and exit"); sp=GetStr(); 
return(UIO_ERROR); 


> 
/* convert float Data to double */ 
for(i=0;1<nDataElems; dDE(ij]=DataElemp[i] , i++); 
rows=cols=0; 


1 heed From the instrument name, extract rows, cols and weight matrix name. 
ia The name of instrument is in the format: 
= "rows(Cany non number character)colsWtname", i.e. 245.12W011 
2 => rows=245, cols=12, Wtname=W011 my 


sscanf(ug_getinstrt(),"%d%cxzd%s", &rows, &c,&cols, inst); 


/* Check for consistency between instrument name and its data size */ 
if (nDataElems! =(rows*cols) )€ 

PutFmt("Instrument nDataPoints/name error. cols=%i, rows=41, mData=%i"', 
rows,cols,nDataElems); 

printf("Instrument nDataPoints/name error. cols=%i1, rows=%1, nData=%i", 
rows ,cols,nDataElems); BEEP; 

error_flag=SET; PutStr("\nPress RETURN to stop and exit"); sp=GetStr(); 
return(UIO_EXIT); 


) 
rewind(mat_i_f); /* to prevent over filling mat_i_f file */ 
/* Here is essentially a small matalb code that calculates and saves the SVDs of 
* the weight matrices. Assuming the name of the weight matrix is WtName, the 


* code is: 
* >> save matif WtName ; 
* >> load matif ; 
* >> sWtName=(sWtName svd(WtName)]; 
* >> clear WtName ; % this is not a must. It is for Wt matrices will not 
* appear when saving the matlab workspace at UIO_Learn_End(). 
* 
/ 
savemat(mat_i_f,1000,inst,rows,cols,0,dDE, (double *)0); /* save */ 


fputs("load ", matlab); fputs(MATLABIF, matlab);fputs(" \n", matlab); /* load */ 

fputs(''s" matlab); fputs(inst,mat lab); fputs("=(s",matlab); 

fputs(inst,matlab); fputs(" svd(",matlab); fputs(inst,matlab); 

fputs("));\n'"',matlab); /* s..=(s.. svd(Wt)]; */ 

fputs("clear ", matlab); fputs(inst, matlab); fputs(" \n", matlab); /* clear */ 
return(ret_val); 

) 


NINT UIO_Rewind() 
{ 


NINT ret_val = UIO OK; 
return(ret_val); 


- 

NINT UIO_ Explain(LayN, nPEs, Datap) 
NINT LayN; 

NINT mPEs; 

SREAL *Oatap; 

C 


NINT ret_val = UIO_ OK; 
return(ret_val); 


232 


[ REAARRERKRAAEREREKKKAEERKEREEEEKKAKEERE RK KREREEKEEKRERKEKKEEKKKERKERREE / 


/* $Author : $ Shahar Dror ae 
‘hel $Date : $ 22 Sep 1991 */ 
| hl $Source >: $ matlab link.h */ 


[RII III ISITE IIHT I IAI II TIS ISI I ISIS SSI AS SSSI ISIN N SO IISA ISS ISIS SS SASOSSS SY 


#include <sys/types.h> 
#include <sys/stat.h> 
#include <stdio.h> 


#define MATLABPIPE "matlab.pipe" 


FILE *matlab; 

FILE *matpipe; 

FILE *mat_i_f; /* matlab interface file */ 

char *MATLABIF="matlabIF"; /* location of interface file for matlab*/ 
char svdn(20J=""s"° /* file name of SVD results*/ 


int start_matlab(¢) 
€ 


time_t *tlc; /* pointer to current time in sec from GMT00:00:00 /’70*/ 
char tstr(5]; /* pointer to hold integer to string conversion */ 
int fd; 
int flags; 
if((matlab = popen("matlab >/dev/null 2>&1", "w'') )==NULL){ 
printf ("error invoking matlab\n");exit (1); 


setbuf(matlab, (char *) NULL); 


jg mknod(MATLABPIPE, S_IFIFO | S_IREAD | S_IWRITE); 
if((matpipe = fopen(MATLABPIPE, "r+''))==NULL){ 
PutFmt("error opening %s",MATLABPIPE); return(UIO_ERROR); 


> 
setbuf(matpipe, (char *) NULL); ul 


sprintf(tstr,"%i", (time(tlc)%100000) ); 
strcat(MATLABIF,tstr);strcat(MATLABIF,".mat"); /* get a unique name */ 
if((mat_i_f=fopen( MATLABIF, "w'') }==NULL )¢ 
PutFmt("error opening %s",MATLABIF); return(UIO_ERROR); 
) 


setbuf(mat_i_f, (char *) NULL); 
return(UIO_ OK); 
) 


void stop_matlab( ) 
£ 


fclose(mat_i_f); 
fputs("delete ", matlab); fputs(MATLABIF, matlab); fputs("\n", matlab); 
fputs("quit\n", matlab); 
pclose(matlab); 
fclose(matpipe); 
unl ink(MATLABPIPE); 
) 


typedef struct ¢ 
long type; /* type */ 
long mrows; /* row dimension */ 
long ncols; /* column dimension */ 
long imagf; /* flag indicating imag part */ 
long namlen; /* name length (including NULL) */ 
> Fmatrix; 


/* Here is an example that uses ‘savemat’ to save two matrices to disk, 


* the second of which is complex: 
* 


ta FILE *fp- 

i double xyz(500], ar(500), ai[(500); 

= fp = fopen("foo.mat", "wb'') >: 

id savemat(fp, 1000, "xyz", 2, 3, 0, xyz, (double *)0); 

“ savemat(fp, 1000, "a", 5, 5, 1, ar, ai); 

~ fclose( fp); 

* See .../matlab/loadsave/savemat.c for more details. = / 
void savemat(fp, type, pname, mrows, ncols, imagf, preal, pimag) 
ELLE. *fo; /* File pointer */ 
int type; /* Type flag: 1000 for Sun */ 
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char “pname; /* pointer to matrix name */ 
int mrows; /* row dimension */ 

int ncols; /* column dimension */ 

int imagf; /* imaginary flag */ 

double *preal; /* pointer to real data */ 
double *pimag; /* pointer to imag data */ 


Fmatrix x; 
int mn; 


x.type = type; 


x.mrows = mrows; 
x.ncols = ncols; 
x.imagf = imagf; 


x.namlen = strlen(pname) + 1; 
mn = X.mrows * x.ncols; 


fwrite(&x, sizeof (Fmatrix), 1, fp); 
fwrite(pname, sizeof(char), (int)x.namlen, fp); 
fwrite(preal, sizeof(double), mn, fp); 
if Cimagf) ¢ 
fwrite(pimag, sizeof(double), mn, fp); 
} 
} 
/* Here is an example that uses ‘loadmat’ to load a matrix from a MAT-file: 
* 
FILES fp; 
char name(20]; 
int type, mrows, ncols, imagf; 
double “xr, *xi; 
fp = fopen("foo.mat","rb""); 
loadmat(fp, &type, name, &mrows, &ncols, &imagf, &xr, &xi); 
fclose(fp); 
free(xr); 
if Cimagf) free(xi); 
* See .../matlab/loadsave/loadmat.c for more details. SF, 
int loadmat(fp, type, pname, mrows, ncols, imagf, preal, pimag) 
ies 


+t 4 &@ 4 He H H 


* 


FILE “fp: File pointer */ 

int *type; /* Type flag: see reference section of guide */ 
int *mrows; /* row dimension */ 

int *ncols; /* column dimension */ 

int *imagf; /* imaginary flag */ 

char *“pname; /* pointer to matrix name */ 


double **preal; /* pointer to real data */ 
double **pimag; /* pointer to imag data */ 


char *malloc(); 
Fmatrix x; 
int mn, namlen; 


if (fread((char *)&x, sizeof(Fmatrix), 1, fp) != 1) ¢ 
return(1); 


} 

*type = X.type; 
*mrowS = X.mrows; 
*ncols = x.ncols; 
*imagf = x.imagf; 
namlen = x.namlen; 


mn = X.mrows * x.ncols; 


if (fread(pname, sizeof(char), namlen, fp) != namlen) return(1); 
if ('(*preal = (double *)malloc(mn*sizeof(double)))) ¢ 
printf("\nError: Variable too big to load\n"); 
return(1); 


di 

if (fread(*preal, sizeof(double), mn, fp) != mn) ¢ 
free(*preal); 
return(1); 


) 


if (x.imagf) ¢€ 
if ('(*pimag = (double *)malloc(mn*sizeof(double)))) ¢ 
printf("\nError: Variable too big to load\n"); 
free(*preal ); 
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return(1); 


} 
if (fread(*pimag, sizeof(double), mn, fp) != mn) {¢ 
free(*pimag); 
free(*preal); 
return(1); 
} 
V 
return(Q); 
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/* file: f18l_modelN.h 

~ date: 8 Feb 1992 

_ revised: 28 March 1992 

id programmer: Shahar Oror 

* 

* This file contains the definitions and the model of the f18 fighter 

* aircraft. Also it defines initial conditions. 

* 

* The names of the functions appearing in this file must remain unchanged! 
* 


‘We te te Se He te te He te te He He te He te te te ee te te te Ke Ke Ke KKK KKK RAEKRRRKRKRKRKURERR AHHH REREEEERERHRRRRERRERUERE / 


/* All definitions in this box MUST appear! ! If don’t care set to default=1(one) 
unless otherwise specified by: "!def=.. = 


/ écansassees\ /seeseesees\ /ssssessess\ /sessessees\ /seesessee=\ /s22s22=22=\ /====*/ 

#define plant_input_order 6 /* order of plant inputs vector */ 

#define plant_output_order 3 /* order of plant outputs */ 

#define in_history 3 /* # of delayed input parameters */ 

#define out_history 3 /* # of delayed output parameters */ 

#define Nmodels 2 /* # of models to train on(>=1). !def=1 

/*®*sssssssses/\sssssrssss/\sssssssssa/\ssssssssss/\ssssssssss/\ssssssss==/\= s=s*/ 

/*a The following definitions are in effect ONLY in this header file. The a 
@ names can be changed, after checking that the desired name does not exista 
@ in other files that are in use together with this header-file !!! a*/ 

#def ine plant_order 25 /* order of plant model =) 


/*qaadaaaaadaadaaadaddadaadadaaaaaadaadaaddaadaadaadadadadadaaadaaadaadaaaaaaaaaa* / 


/* Declarations */ 

/*® ssssssssssssss */ 
/* The values under declarations MUST appear!! If "don’t care" set to default, 
which is 0 (zero) in the correct type, unless otherwise specified by "def=..."*/ 
/*sszsssssssa\/ssasassssa\/sssssssasa\/ssssesssas\/stssssssas\ /sssssssssa\/ss==*/ 


double ts=0.02; /* Sampling interval for discretization. def=1.0 */ 
double tau=17.74; /* Characteristic(longest) time constant [sec] (Phugoid)*/ 

/* Frequency band of interest [Hz] in the current model */ 

double Fmax=5.0; /* Upper frequency (used for swept wave test) */ 

double Fmin=0.005; /* Lowest frequency (used for swept wave test) */ 

double excite (plant input_order]=(0.0);/* Plant input commands = NN outputs */ 
double bias_I (Nmodels] (plant_input_ order] =(0. 0); /*Plant input bias matrix*/ 
/¥aessseessa/\ssescessss/\cetssseses /\ssesseeess /\seessesees/\sssesseess/\s===*/ 


/*a The following declarations are in effect ONLY in this header file. The @) 
@ names can be changed, after checking that the desired name does not exista 
@ in other files that are in use together with this header-file !!! a*/ 


double excite ee Sede t ante input_order] = /* Scaling the input commands */ 
(10,1.0,1.0, 1.0, 1.0.1.0): 

/* Initial conditions. If !=0 put ‘them here! ine 

double plant_state_ic(plant_order]=(0.0);/* Plant initial conditions vector */ 


double plant_A(Nmodels] (plant_order*plant_order];/* Plant dynamic matrix */ 
double plant_8 (Nmodels] (plant_order*plant_input_order];/* Plant control matrix*/ 
double plant_C(Nmodels] (plant_output_order*plant_order];/*Plant measure matrix*/ 
double plant _0 (Nmodels) (plant_ output_ _order*plant_ input_ order]; /*Plant ff matrx*/ 
double bias_O(Nmodels] (plant_output_order]=(0.0);/*Plant output bias matrix*/ 


/*dddddAaAadAdAdAAdAdAAAAdAAdAAIAdAddAdAdAddAdAdAAddddAdAddAAdAddAdAddddAdAdAadAaAadAdAadAdAAAA* / 


load_model (mnr ) /* Loading matrices of current model to work with */ 
int mnr; /* Model number, 1,2,...,Nmodels */ 

{ 

double *xr,*xi; /* Real and imaginary parts of the loaded matrix */ 


int type, mrows, ncols, imagf; 
int ret_val=UIO_OK; 


ahem legis /* Indexes */ 
char name (20); /* Name of loaded matrix. Not used. */ 
char intstr(10]; /* String to hold integer to string conversion */ 


char matfile(25]); /* string to .mat file name for model matrices */ 
FILE *ck_open_file(),*fp; 
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sprintfCintstr,"%i",mnr); 
streat(intstr,".mat")> /* Create extension to model matrices as XXX#.mat */ 


strepy(matfile,"ZzFld");strcat(matfile,intstr); 
if((fp=ck_open_file(matfile,"r'))==NULL) return(UIO_ERROR); 
loadmat(fp,&type,name, &mrows, &ncols ,ximagf,&xr,&xi); fclose(fp); 
if((mrows != plant_ order) ir ‘(ncols {= plant_ ‘order)){ 
PutFmt("\n!!Plant mismatch. plant_A (model ANE SGT eA plant_order=4i", 
mnr,mrows,ncols,plant_order); 
PutStr("\nError in load_model(). Press RETURN to stop and exit"); 
GetStr(); return(UIO_ERROR); 
} 


for(i=0;i<mrows*ncols;plant_A(mnr-1) Ci) =xr (Ci), i++); 


strepy(matfile,"zGld");strcat(matfile, intstr); 
if((fp=ck_open_file(matfile,"r"))==NULL) return(UIO_ERROR); 

loadmat( fp,&type, name, &mrows ,&ncols,&imagf,&xr,&xi); fclose(fp); 
if((mrows != plant_ order) a ‘(ncols i= plant_ input_order)){ 
PutFmt("\n!Plant mismatch. plant_B(model%i]is(%ix%T) ,but set to (%ix%i]" 
»mnr,mrows,ncols,plant_order »plant_ input_order); 
PutStr("\nError in load_model(). Press RETURN to stop and exit"); 
GetStr(); return(UIO_ERROR); 


for(i=0;i<mrows*ncols;plant_B{mnr-1) CiJ=xr (i), i++); 


strepy(matfile,"ZzHld");:strcat(matfile,intstr); 
if(( fp=ck_open_file(matfile,"r'))==NULL) return(UIO_ERROR); 

loadmat( fp,&type,name, &mrows ,&ncols,&imagf,&xr,&xi); fclose(fp); 
if((mrows != plant_output_order) ;{| (ncols != plant_order))¢ 
PutFmt("\n! Plant mismatch.plant_C(model %i)is(%ix%i],but set to (%ix%i]" 
»mnr,Mrows ,ncols,plant_output_order,plant_order); 
PutStr("\nError in load | model (). Press RETURN to stop and exit"); 
Getstro return(UIO_ERROR); 


} 
for(i=0;i<mrows*ncols;plant_C(mnr-1) CiJ=xr(li], i++); 


strepy(matfile,"zDld"):strcat(matfile,intstr); 
if((fp=ck_open_ file(matfile, "r))==NULL) return(UIO_ERROR); 
loadmat( fp, &type, name, &mrows, &ncols ,&imagt, &xr,&xi); fclose( fp); 
if((mrows != plant _output_ order) - (ncols 1s plant_input_order))¢€ 
PutFmt("\n!Plant mismatch.plant_D{model %i]is(%ix%i],but set to (%ix%i]" 
»mnr,mrows,ncols,plant_output_order ,plant_input_ order); 
PutStr("\nerror in load | model (). Press RETURN to stop and exit"); 
GetStr(); return(UIO_ERROR); 


for(i=0;1<mrows*ncols;plant_D(mnr-1) Ci] =xr Ci], i++); 


strepy(matfile,"zBil");strcat(matfile,intstr); 
if ((fp=ck_open_ file(matfile, Wp )s =NULL)¢€ 
PutStr("\nAssuming no input -biasing*G"); 
printf ("\nAssuming no input-biasing*G"); 
> 


else( 

loadmat( fp, &type, name, &mrows, &ncols imagf, &xr,&xi); fclose(fp); 
if! (C(mrows= =plant_ input_ order)! “(ncols == plant_input_order)) 

&& ((mrows == 1)}!(ncols == 1))))¢€ 
PutFmt("\n!Bias mismatch. bias_ I (model “%i]is(%ix%i],but set to (%ix%i)" 

,mnr,mrows ,ncols,1,plant_input_order); 
PutStr("\nError in load _model(). Press RETURN to stop & exit"); 
GetStr(); return(UIO_ERROR); 
d/*end if*/ 
for(1=0;i<mrows*ncols;bias_I (mnr-1] Ci]=xrCi),i++); 

)/*end else*/ 


strepy(matfile,"ZzBol");strcat(matfile,intstr); 
if (( fp=ck_open_ file(matfile, "P1) )==NULL) € 
PutStr("\nAssuming no output-biasing*G"); 
printf("\nAssuming no output-biasing*G"); 
) 
else{ 
loadmat( fp, &type, name, &mrows, &ncols Gimagt xr, &x1); fclose( fp); 
if¢!(((mrows==plant_ output_ order)! i(ncols ==plant_output_order)) 
&& ((mrows == 1)!!(ncols == 'N))))€ 
PutFmt("\n!Bias mismatch.bias_Ol(model %i]is(%ix%i],but set to (%ix%i)" 
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,mnr,mrows,ncols,1,plant_output_order); 
PutStr("\nError in load_model(). Press RETURN to stop & exit"); 
GetStr(); return(UIO_ ERROR); 
d/*end if*/ 
for(i=0; i<mrows*ncols;bias_O(mnr-1] Ci]=xr(i], i++); 
}/*end else*/ 


return(ret_val); 


PLANTCic,mnr,Plant_input,Plant_out)/*Mathematical description of plant output */ 


int ic; /* (Initial) conditions type indicator. IClearn,ICtest, TEST or LEARN */ 
int mor; /* Model # in effect for current desired output calculation */ 

/* 1 - non damaged, original. >=2 - impaired, other */ 
double Plant_input (plant_input_order]; /* Control inputs vector */ 


double Plant_out(plant_output_order]; /*f18 measurements+model coding */ 


f 
/* This function is written in details in order to be able to put in any non-linearities at any 
point in the calculation*/ 


static int mnr_last=1; /* previos model number. Initially the first. */ 
static double plant_state(Nmodels] (plant_order]; /* £18 state-vectors */ 
double plant_msr_state(plant_output_order]; /* £18 msr vector */ 
double wv(plant_order]; /*working vector*/ 

int 1,J; 


I he! First X(k+1) is calculated using X(k), UCk), then - Y(k+1) using X(k+1). */ 
/* For initial condi tions don’t update the states */ 
if(ic==IClearn |, ic==ICtest)¢ 
for (i=0; i<plant_order; i++) 
plant_ state(mnr-1] Li]=( Cic==ICtest) ? 0.0 : plant_state_ic(i]); 
goto [Cjump; 
> 


/* For continuous RBS input (=1) start new model with last state value of 
previous model. */ 
ifCic==TEST1 && mnor!=mnr_last)¢ 
for(i=0; i<plant_order; i++) 
plant_state(mnr-1] (iJ=plant_state(mnr_last-1] Ci]; 
mor_last=mnr; 


/* Plant state update*/ 
for(1=0;1<plant_order; i++){ 


wv(1]=0.0; 
for(j= 0; j<plant_order; j++) /* Ax */ 
wv(i]+= (plant_ A(mnr- 1] Litplant_ order* j]*plant_state(mnr-1] (j]); 
for(j=0; ;j<plant_ input_order; j++) /*+BCutbi) */ 


wVLi]+=(plant_ B (mnr - 17 Ci+plant_ order* j] 
*(Plant_input(j]+bias_I (mnr-1] (j])*excite_amplitude(j]); 


} 
for(1=0; i<plant_order;plant_state(mnr-1] Ci]=wv(i] , i++); LOX Ce) *7 
IC jump: 


/* Plant measurement update*/ 
for(i=0; i<plant_output_order; i++){ /* Cx+Dutbo */ 
Plant_out (i]=0.0; 
for ( j=0; j<plant_order ; j++) 
Plant_out(i]+= 
(plant_C(mnr-1] Citplant_ output_order* j]*plant_state(mnr-1] (j]); 
/* in my case = 0 for( j=0; j<plant_ input_order; j++) 
Plant_out(i]+=(plant_ 'D (mnr - 1] Ci+plant _output_ order* j] 
*(Plant input(j]+bias_ I (mnr-1] Cj] )*excite_amplitudelj]);*/ 
Plant_out(i]+=bias O(mnr-1] (1); 
> 
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eg File: net_paramNtp.h ah 
/* Author: Shahar Dror ay 
/* Date: 9 Mar 92 af 
j* Revised:17 Jun 92 ia 
[es Subject: This header file contains declarations and */ 
lig definitions of parameters that are used in the USRIO. ie A 
ihe It comes before net_util.h and after ref_modelN.h, */ 

[* and uses parameters/values declared there. olf 
This header file also contains functions*/ 

/* used by the USRIO program for its execution. Parameters */ 
/* declared here are global, and accessible to every tld 
/* function in the USRIO, and also to every function a J 
/* hereafter. Declarations and definitions here should not */ 
/* be deleted. ial 
j* This file creates the base parameters that describe the */ 
Ves problem being worked on. Their VALUES may be altered ‘edi | 
/* (but not their names). Default values (i.e if "don’t- *=/ 
/* care") are zero (0), unless otherwise specified. 3/6/92 */ 
ai * 


[ RARER RRRREREEREREREEEEEEEREREREEEREEUEEEEEREE WHHKEKRKRHRREREKRKKEEKREERERKEEE / 


/* Essential definitions. Provide work base for the current problem. !Default=0* 


* The variables are usually self-explanatory by their names. =) 
pene aa aa aa SaaS aaa eee e one e aaa aaa ee ee ae Sea eeSaaSSsssssssssssssssssss=es*/ 
#if Task==classification /* Set up classification mode */ 
#define forced_excitation_length plant_input_order 
#define unforced_excitation_length plant_output_order 
#define type_length 0 
#define excitation_length (forced_excitation _length+unforced_ excitation_length) 
#def ine regressed_ Tnput_ Length excitation _tength 
#def ine input_lay_ Length regressed_ Tnput_ length 
#def ine regressed_input_history_length regressed_input_length*in_history 
#def ine dynamic_output_length 0 
#define regressed_output_length 0 
#define static_output_length Nmodels 
#def ine output_lay_length Nmodels 
#define regressed | output_history_length regressed _output_length*out_history 
#define effector_ltay_length 0 
#endi f 
#if Task==emulation /* Set up emulation mode */ 
#define forced_excitation_length plant_input_order 
#define unforced_excitation_length 0 
#define type_tength Nmodels 
#def ine excitation_length forced_excitation_length 
#def ine regressed_ Tnput_ Length forced_excitation_length 
#define input_lay_length regressed_ input_ length+type_ length 
#def ine regressed_ input_history_tength regressed input_length*in _history 
#def ine dynamic_output_Tength plant output order 
#define static_output_Tength 0 - ms 
#define regressed_output_length plant_output_order 
#define output_ lay_ Length plant_output_order 
#def ine regressed_output_history length regressed_output_length*out_history 
#define effector_lay_tength 0 
#endi f 
#if Task==combined_cl_em /* Set up classification and emulation mode */ 
#define forced_excitation_length plant_input_order 
#define unforced excitation _length 0 
#define type_ Length Nmodels 
#define excitation _tength forced excitation length 
#define regressed _Tnput_ltength forced_excitation_length 
#define input_lay length regressed_ input_ Length+type_ length 
#define regressed_input_history_length regressed_ Input_length*in_ history 
#define dynamic_output_length plant_output order 
#define static_output_Tength Nmodels ~ 
#define regressed_output_length plant_output_order 
#define output_lay_length plant_output_order+type_length 
#define regressed output _history_length regressed output length*out history 
#define effector_Tay_length 0 4 ee - 
#endi f 
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#if Task==controlling /* Set up controlling mode */ 
#define model_input_length 6 

#def ine model_output_length plant_output_order 
#def ine forced_ excitation_length 


plant_input_order 
#def ine unforced_ excitation_length 0 
#define type_ length (model_input_length+model_output_length+Nmodels) 
#define excitation_length forced_excitation_length 
#def ine regressed_input_length excitation_length— 
#define input_lay_length regressed _input_length+type_length 
#define regressed_input_history_length regressed_input_length*in_history 
#define dynamic_output_Tength plant_output_order 

0 


#define static_output_length 

#define regressed _output_length plant_output_order 

#define output_lay_length plant_output_order 

#define regressed output_history_length regressed | output_length*out_history 


#define effector_lay_length model_input_length 


#endi f 


#if Task==cntrlThruPlant 
#define ref_in_length 
#define ref_out_length/*eq plant_out*/ 
#defi.ie plant_in_length/* j; “i 
#define plant_out_length/*eq ref_out*/ 
#define forced_ excitation. length 
#define unforced excitation _ length 


/* Set up controlling mode */ 
ref_input_order 
ref. _output_order 
plant_input_order 
plant_output_order 
ref_input_order 
ref_output_order 


#define type_ length (plant_input_order+plant_output_order ) 
#define excitation_length (forced_excitation _length+unforced_ excitation_length) 
#def ine regressed_ Tnput_ length “(excitation_ length+type_ length) 
#define input_ lay_ length (excitation_ length+type_length) 

#def ine regressed_ input_history_length regressed_input_ length*in history 
#define dynamic_output_length plant_input_order 

#define static_output_ Tength 0 

#define regressed |_output_length 0 

#define output_lay_length dynamic_output_length 


#define 
#define 
#endi f 


regressed_ output_length*out_history 
output_lay_length 


regressed output_ history_length 
effector Vay_ length 


/ bss If the regressed values are not consequent and as described here,* 
* indicate the sorting in the functions regress _input_sort(), and bd 
* regress_output_sort(). These appear down in this header file. is! 


= The user must take care and connect the relevant input nodes to * 

* the command regression layer in the network (while building it)!! */ 
/P®SSSSSSSSSSRSSAAAAAASSSAASSSSSRSSAAASSSSSSSSSSSSS SSS SSS SSS SS SSSSSSSsesSSs=s===*/ 
#define RibsUd 4 /* RBS update divisor */ 

#define NErrConv 20 /* # of converged RMS error for best network check */ 
#define NetSaveMode 0 /* O=binary, 1=ascii, 2=annotated, 3=screen dump */ 
#define ninstElmnt 8190 /* max # of elements sent by an instrument */ 

#define maxInstName 25 /* Max length of instrument name */ 

#define Ninput_types 6 /* # of input types (see *input_name -1) */ 

#de fine maxNsweep_ change 700 /* Max # of sign switching for H->L swept wave*/ 
double dpl (]=(0.5,1.5,2.5,5.0); /* dipole switching times for 0,+1,-1,0,end. */ 
double convergence=2e- => /* Total error RMS convergence criterion */ 


double h=0.01; 
float min_lim=-2.0,max_lim=2.0; 
char NNrslt(20]="NNout.dat"; 


Spatial step size for Jacobian calculation */ 
Limit values for network outputs */ 
file name of NN results in TEST mode*/ 


char DesO(20]="DesOut .dat"; /* file name of desired outputs in TEST mode*/ 
char Inp(20]="Input.dat": /* file name of inputs in TEST mode*/ 
char *input_name[(]=¢ /* '!i!tcheck (#define Ninput_types) up here!!! */ 
“Illegal Input", /* input=0 xy 
“Random Binary", /* input=1 */ 
"Sequential Random Binary", /* input=2 */ 
“Swept H->L Square Wave", /* input=3 */ 
"Swept L->H Square Wave", /* input=4 */ 
"Swept L->H Sine Wave", /* input=5 */ 
"1 sec dipole"); /* input=6 */ 


int input=1; /* Input type indicator, initially set to RBS (=1) */ 
int cmRC=1000; /* Divisor of RBS random TEST cycles model switching */ 
int maxSW=2,minSW=1; /* Divisor,bias of RBS random cycles model switching */ 
int nRTC,ResetTC=10000; /* Reset-training-cycles mode divisor */ 
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/* Place here global parameters that are needed for the specific implementation. 
Check that they don’t exist in the program, including all include files. 
These parameters are to be used within functions present in this 
header file (net_paramXxxX.h). “Wy / 

[Boer creer eee Vo-cereeeree- Vnnnrreenen-- Ye V------------ V----------- */ 

#if Task==cntrlThruPlant 

double refout (Nref_models] [ref_out_length] ; 

double plantout [Nmodels] [plant_out_length]; 

double outerr([plant_out_length]; 


double W(plant_out _tength}= ={(1.0,1.0,1.0); /* Output weighting vector */ 
#endi f 

j* Oo Oe eee eee a ee a Woe eannaewaaae# CoN Len Yen a ee NeoGeo Bbloleeee = = 
regress _ input_ sort(RegInSort) /* Input sorting function */ 

int RegInSort [regressed_ input_length] ; 

{ 

int i; /* Index counter */ 


/* Default sorting is in sequential order */ 
for(i=0;i<regressed_input_length;RegInSort [i] =1, i++); 


regress_output_sort(RegOutSort) /* Output sorting function */ 
int RegOutSort (regressed_output_length+1]; 
¢ 


mnt. 13 /* Index counter */ 

/* Default sorting is in sequential order */ 
for(i=0;i<regressed_output_length;RegOutSort [i}=i1,1++); 

) 


int Limits(mnr, out_val) /* Returns 1 if inside limits. Else returns 0 */ 
/* User can write any limit function, per each output, or common to all outputs, 
with min/max or abs values, etc. */ 
int mor; /* Model # in effect for current Jacobian calculation 
1 - non damaged, original. >=2 - other */ 
float out_val [output_lay_length]; 
{ 
int i; /* counter */ 
int lim; /* if=1, within limits, if=0, out of limits */ 
for(1=0;(i<output_lay_length && lLim!=0); i++) 
Lim=(Cout_val (iJ>=min_lim && out_val [i] <=max_lim)?1:0); 
return(lim); 
> 


jacobian(mnr ,plant_input,plant_output, jacb) /* Jacobian calculation */ 
double jacb[plant_output _order*plant_ input_order]; /* Jacobian matrix */ 
double plant_input [plant_input_order]; /* Control inputs vector */ 
double plant_ _output [plant _output_ order]; /* Current plant outputs vector */ 
Int mor; /* Model # in effect for current Jacobian calculation 

1 - non damaged, original. >=2 - other */ 


{ 
static double old_plant_input [Nmodels] [plant_in_length}={0.0);/* Last inputs */ 


double dh; /* Spatial increment including sign */ 
double plant_ out (plant_output_order]; /* plant output @ perturbated inputs */ 
int i,j; /* indexes */ 


for(i=0;1<plant_input_order; i++){ 
Poi Choose the sign of | increment to ensure staying in the domain */ 
dh=copysign(h, (old_plant_input (mnr-1] Ci] -plant_input li} )); 
plant_ input [i]+=dh3 
PLANT (JACOB, mnr ,plant_input,plant_out); 
for(j=0; j<plant. _output_ order; j++) 
jacb[i*plant-output_order+j}= (plant_out(j]-plant_output [Jj] )/dh; 


} 
old_plant_input (mnr-1] (i]=(plant_input [i] -=dh); 
} 


> 

desired_output(ic,mnr,Plant_input,dso) /* Get desired output */ 

int 1c? /* (Initial) conditions type indicator. IClearn, ICtest, TEST, LEARN, 
TEST1 or JACOB to indicate Jacobian calculation. */ 

Int mnr; /* Model # in effect for current desired output calculation */ 
/* 1 - non damaged, original. >=2 - impaired, other */ 

double Plant_input [forced_excitation_length]; /* Control inputs vector */ 


double dso[output_lay_length]; 
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{ 
extern double Effector [(Nmodels] (effector_lay_length+t1]; 


double jac(plant_output_order*plant_input_order]; /* Jacobian matrix */ 
extern int error_flag; 
int 4, 12 /* Index counter */ 

if (Task==emul ation) PLANTCic,mnr,Plant_input ,dso); 


else if (Task==combined_cl_em){ 
PLANTCic,mnr,Plant_input,dso); 
for (i=0; i<type_length; i++) 
dso (dynamic_output_length+i] =(i==(mnr-1)?1.0:-1.0); 
BCD(mnr , &dso[dynamic_output_length] ); 


else if(Task==controlling)¢€ 
REFCic,mnr,Plant_input,dso); 
PLANTCic,mnr, &Effector (mnr-1] (0] ,&plantout (mnr-1] [0] ); 


) 
else if(Task==cntrlThruP lant )¢ 
/* In this implementation, for every plant there is a corresponding 
reference. If the situation is not so, there should be a matching 
between the plant # (mnr) and the value sent to REF(). */ 
REFCic,mnr,Plant_input, &refout (mnr-1] [0] ); 
PLANT(ic,mnr,&€ffector (mnr-1] [0] ,&plantout (mnr-1] [0] ); 
for(i=0;i<plant_out_length; 1++) 
outerr (1]=refout (mnr-1] Ci] -plantout (mnr-1] [i]; 


jJacobian(mnr , &Effector (mnr-1] (0] ,&plantout (mnr-1] (0), jac); 
for(i=0; i<plant_in_length; i++)¢ 
dso [iJ]=0.0; 
for ( j=0; j<plant_out_length; j++) 
dsoli]+=(W[j]*outerr(j]*jac(i*plant_out_length+j]); 
dso(i]+=Effector (mnr-1] [i]; 


) 
) 
/* User put before this line more tasks if desired. The “else must be last! */ 
else¢ 


PutFmt("\n!!No desired_output() for specified Task (= %1)", Task); 
printf("\n!!No desired _output() for specified Task (= %i1)",Task); 


error_flag=SET; BEEP; 

) 
) 
unforced_excitation(ic,mnr,Plant_input,ue) /* Get unforced excitation */ 
int 1c; /* (Initial) conditions type indicator. [Clearn, TEST or LEARN */ 
int mnr; /* Model # in effect for current desired output calculation */ 

/* 1 - non damaged, original. >=2 - impaired, other */ 

double Plant_input(forced_excitation_length] ; /* Control inputs vector */ 
double ue[unforced_excitation_length*1]; 
{ 


static double ueOld(Nmodels] (unforced_excitation_length+1]; 
extern double desoutOld(Nmodels] [output_lay_length]; 
extern int error_flag; 
int i; 
if (Task==classification){€ 
for(i=0;i<unforced excitation length;ue[i]=ueOld(mnr-1] Ci], i++); 


r 


PLANT(ic,mnr,Plant_input ,&ueOld(mnr-1] (01); 


) 
else if(Task==controlling){/* The prior value of the reference model */ 
/*for(i=0; i<unforced_excitation_length; i++) 
ue [i] =desoutOld(mnr-T] [i] ;*/ 
, 
else if(Task==emulation); /* No unforced_excitation in emulation*/ 
else if(Task==cntrlThruP lant )¢ 


for(i=0; i<unforced_excitation_length; i++)¢ 


ted ifCic==IClearn |! ic==LEARN)*/ 
ue(iJ=refout(mnr-1]) [Ci]; /* it is the old value */ 
Ves ifCic==ICtest |} ic==TEST) 


ue{iJ=plantout (mnr-1] (i];/* it is the old value */ 
) 
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/* User put before this line more tasks if desired. The "else" must be last! */ 


else{ 
PutFmt("\n!!No unforced_excitation() for specified Task (= %1)",Task); 


printf("\n!!No unforced _excitation() for specifted Task (= %1)",Task); 


error_flag=SET; BEEP: 
} 
) 
Type(ic,mnr,nnout,effect_out ,mtype) /* Model type identification */ 
lint ice /* Cinitial) conditions type indicator. IClearn, TEST or LEARN */ 
int mnr; /* The model # */ 


double nnout foutput_lay_length]; /* Network output. Not always needed. */ 
double effect_outleffector_lay_length+1];/* Effector value. Not always needed.*/ 
double mtype[type_length+1]; /* Its representation */ 

{ 


extern int error_flag; 
int i; /* Index counter */ 
if (Task==emulat ion) 
for(i=0;1<type_length;mtype(i]=(1==(mnr-1)71.0:-1.0), i++); 
else if(Task==controlling)¢€ 
for(i=0;i<effector_lay_length;mtype(li]=effect_out [i], i++); 
for(i=0; i<output_lay_length; i++) 
mtype {iteffector_lay_length] =nnout [i]; 
for(1=0;1<Nmodels; i++) 
mtype[i+type_length-Nmodels]=(1==(mnr-1)71.0:-1.0); 


else if(Task==cntrlThruPlant)¢ 


for(i=0;i<effector_lay_length; i++)¢ 
mtype [il=effect_out[i]; 


} 
for(i=0;1<plant_out_length; i++) /* it is the old value */ 
mtypeCi+effector_lay_length]=plantout (mnr-1] [i]; 


} 


/* User put before this line more tasks if desired. The "else" must be last! */ 


else{ 
PutFmt("\n!!No Type() for specified Task (= %i)",Task); 
printf("\n!!No Type() for specified Task (= %i1)",Task); 


error_flag=SET; BEEP; 
} 
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[MRARRRRARERER RERREREEERARNE RENAE EES ee ee 


7" File: net_util.h as, 
/* Author: Shahar Dror ay 
V ha Date: 9 Mar 92 */ 
Thad Revised: «/ 
/* Subject: This header file contains declarations and aid 
/* definitions of parameters that are used in the USRIO. ial 
Vas It comes after net_param.h, and uses parameters/values Kl 
/* declared there. This header file also contains functions*/ 

jie used by the USRIO program for its execution. Parameters */ 
Fabs! declared here are global, and accessible to every */ 
/* function in the USRIO, and also to every function */ 
j* hereafter. Declarations and definitions here should not */ 

isi be deleted. ial 
te ys 


[RARER RRARARARARAENE RARE ERNE eR ee Ne ee ee ees 


#define RBS (2.0*(random( )&01)-1.0) /* Random Binary Sequence */ 
#define DIPOLE(T) (T>=dpl (0]&&T<dpl (1) )?1.0:¢¢T>=dpl (1] &&T<dpl (2] )?-1.0:0.0) 
#define Large number Cint)(Cexp2(31.0)-1) /* Large integer # */ 
#define PI 3.174159265358979 /* The constant pi */ 

double desout (Nmodels] (output_lay_length]; /* Desired output */ 

double desoutOld{(Nmodels) (output _lay_ length] ; /* Last desired output */ 


double unfexc(Nmodels] (unforced excitation _length+1]=(0.0);/*Unforced excita’n*/ 
double ComReg (Nmodels] [regressed_ input _ history_ length] ;/*Command regres output*/ 
double OutReg (Nmodels] (regressed output_history_ length+1]; /*Out-regress output*/ 
double Effector (Nmodels] (effector _lay_length+1]=(0.0); /*Effector/Controler out*/ 
double input_data(Nmodels] Cinput_ Tay_lTength] =(0.0); /* Input-layer data */ 


double err_RMS,max_RMS=0.0,min_RMS=1. 203 /* Relative mean square output error */ 
double dmin, df, f. 1: /* Log10(Fmin), deltaF,F and Time for swept wave */ 

double df log, df lin, fflog, flog, flin; /* Frequences for swept wave */ 
double one=1.0; /* unit value for H->L swept square wave */ 

double min_sum, sum(NErrConv+1]={0.0); /* RMS error convergence parameters */ 


char *BestNetName; /* file name of the best network results lf 

int effector_lay= =7; /* Controller output layer #4 (initial guess) */ 

int rand_inputs (forced excitation length] ;/*Array of random unrepeated inputs */ 
int rand_models (Nmodels+1]; /* Array of randomly unrepeated model #s. */ 
int RIsort (regressed_input_length]; /*Ordering input units for regression*/ 
int ROsort (regressed_output_length+1]; /*Ordering output units for regression*/ 
int Tau; /* the value "tau" measured in ts units ==tau/ts*/ 

int nlayp,ninp,noutp; /* Pointers to: # of layers,# of input/output nodes */ 
int history; /* # of cycles w/o learning to fill regression layers */ 

int inTest; /* # of input under test, or training input in LEARN */ 

int model_nbr; /* # of model under test, or training */ 

int mRC; /* minimum Recall/test/learn-switching Cycles */ 

int MinErTc ,MaxErTC; /* min/max error training/testing cycles indicator */ 
int conv_test (2]= =(0,1); /* Convergence test counters */ 

int Nsweep_change,sweep_count; /* High->Low swept square wave counters */ 

int sweep_change_ times [(maxNsweep | change]; /* H-L switching times vector */ 

int netSave _flag= RESET; /* Flag when (RESET)/SET (don’t)/save best network */ 
int TestDataSave_flag= =RESET:; /* Flag (CRESET)/SET (don’t)/save TEST results */ 
Int error_flag=RESET; /* Flag when (RESET)/SET an error (didn’t) occurred */ 
int intrrpt_flag= SET: /* Flag when (RESET)/SET run was (not) interrupted */ 
int Almodel flag= SET: /* Flag when (RESET)/SET (more than) one model active */ 
int Rstart _flag= RESET; /* Flag when (RESET)/SET (not) passed UIO_Recall_ Start*/ 
int ResetTrain _flag=SET; /* Flag when (RESET)/SET (dont) use Reset-trainig mode*/ 
int ResetLimit_ flag= RESET; /*(RESET)/SET (dont) Reset when output exceed Limits*/ 
int netReset _flag= RESET; /*CRESET)/SET (dont) Reset network */ 

int eflay_ f lag=RESET; /* Flag when (RESET)/SET (dont) recheck effector_lay #*/ 
FILE *NNo;  /* NN output file in TEST mode */ 

FICE *DESo; /* Desired output file in TEST mode */ 

FILE *InPut; /* Input data file in TEST mode */ 


/* In order to prevent compiling errors */ 
#ifndef Nref_models 

#define Nref_models 0 

load ref _model(){} 

#endi f 
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int check_net() 
/* This routine checks that the # of input/ouput PEs match the desired in the USRIO program. If 
there 1s a match the routine returns UIO_OK, on failure it returns UIO_ERROR */ 


{ 
/* Network parameters */ 
int ltype; /* Network parameters */ 


static int nlays=0; /* # of layers reminder */ 

char *csp,*netnp; /* Network pointers*/ 

char “sp; /* pointer to a string */ 

long ctrp(8]; /* An array of 8 to interact with control-strategy*/ 


/* Read network parameters */ 
ug_rdnetinf( &nlayp, &ninp, &noutp, &ltype, &csp, &netnp ); 
/* Check if the built network fits the current plant, control and history */ 
if(ninp != input_lay_length)¢ 
PutFmt('"\n!! The network doesn’t match. nPEs_in=%i, but has to be :%i", 
ninp, input_lay_length); 
PutStr("\nPress RETURN to stop and exit"); GetStr(); 
return(UIO_ ERROR); 
iy 
if(noutp != output_lay_length){ 
PutFmt("\n!! The network doesn’t match. nPEs_out=4%i, but has to be :4i", 
noutp, output_lay_length); 
PutStr("\nPress RETURN to stop and exit"); GetStr(); 
return(UIO_ERROR); 
) 
/* Check consistency between parameters and network structure */ 
if(nlayp<5 && regressed_output_length>0)¢ 
PutFmt("\n!! OutReg layer requested. At least 5 layers needed. nLays=Ai" 
snlayp); BEEP; 
PutStr("\nPress RETURN to stop and exit"); GetStr(); 
return(UIO_ERROR); 
y 
/* Set Effector output layer #, and specify it to the control strategy via AUX2 
counter */ 
if(effector_lay_length>0 && (nlays!=nlayp |} eflay_flag==SET))¢ 
PutFmt("\nEnter effector-layer # (input_lay=0) (%i] ",effector_lay); 
sp=GetStr(); 
if(*sp) sscanf(sp,"%i",&effector_lay); 
ug_ctrio(ctrp,0); /* read counters content */ 
if(effector_lay==(nlayp-1)) ctrp[4]=0; 
else ctrp[4]= (long)effector_lay;/* Set Aux2= effector- layer # */ 
ug ctrio(ctrp,1); _/* write counters content */ 


nlays=nl ayp; /* Set #4 of layers reminder to # of layers */ 
efLay_flag==RESET; /* Reset effector_lay # checking */ 
y 


/* save network name for storing best network results */ 
BestNetName=strcat(netnp,"Bst"); 


return(UIO_ OK); 
} 


void set_net() 

{ 

/* Setting and initializing different network and other related parameters */ 
/* THIS FUNCTION MAY BE EDITED (BY ADDING!!) BY USER */ 


Tau=(int)(ceil(tau/ts)); /* Setting time constant-tau in units of ts */ 

/* Sort the # if the output units to be regressed if not consecutive */ 
regress output _sort(ROsort); 

/* Sort the # if the input units to be regressed if not consecutive */ 
regress_input_sort(Rl sort); 
) 


FILE *ck_open_file(fpn, type) 
char *fpn, *type; 
fal Opens a file for r/w . 
fpn - is the name of the file to be open. 
type - is "r" for read, "w"' for write. 
On success, it return the file pointer. If the file can not be open it 
gives a statement and returns NULL. ail 
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FILE *fp; 
1f(( fp=fopen( fpn, type) )==NULL ){ 
printf('"' Can’t open (for %s) file: %s", type, fpn); 
PutFmt('"* Can’t open (for %s) file: %s", type, fpn); 
) 

return( fp); 

) 


int open_TEST_files() /* Open files for saving TEST results */ 
char *sp; 


PutFmt("\nEnter file name to save Inputs [%s] ",Inp); 
sp=GetStr(); 

if(*sp) sscanf(sp,'*%s", Inp); 

if((InPut=ck_open_file(Inp,"w'!) )}==NULL) return(UIO_ERROR); 
PutFmt("\nEnter file name to save NN output [%s) ",NNrslt); 
sp=GetStr(); 

if(*sp) sscanf(sp,'*%s",NNrslt); 

if((NNo=ck_open_file(NNrslt,"w''))==NULL) return(UIO_ERROR); 
PutFmt('"\nEnter file name to save desired output (%s] ",DesO); 
sp=GetStr(); 

if(*sp) sscanf(sp,'"%s",DesO); 

if ((DESo=ck_open_file(DesO, "w"") )==NULL) return(UIO_ERROR); 
return(UIO_ OK); 
} 


void close_files() 


Vid Closes all opened files bei f 

{ 

fclose(InPut); /* Closing input file Cif exists) */ 
fclose(NNo); /* Closing NN output file (if exists) */ 
fclose(DESo); /* Closing desired output file (if exists) */ 
return; 

) 


void rand_array(N,rnd_vec) 
/* This function generates unrepeated (N) integer random #s valued 1 to N. */ 


int N> /* # of elements in the random vector.(==forced_excitation_length)*/ 
int *rnd_vec; /* Array containing unrepeated random #s */ 

{ 

lite Wale 


rnd_vec (0) =(random( )2N )+1; 
for( j=1; J<N; j++) 
do{ /* generate new non repeatetive random input */ 
rnd_vec [j]=(random( )2N )+1; 
for(i=0;¢(i<j && rnd_vec(j) != rnd_vec(i));1++); 
} whileci != J); 


y 
y 
BCD(N,Nbcd) /* Convert N to its Binary Coded Decimal representation */ 
Int N; /* The # to be converted */ 
double Nbcd[type_length+1); /* Its BCD representation */ 
{ 
Int 1; /* index */ 
int ne; /* copy of N */ 


nc=2*(N-1); 
for(1=0;i<type_length;Nbcd[i)=(double)( (nc=nc/2)%2)*2-1, i++); 
7 


/* From here on are function of different input-commands for testing. */ 
y* we 


int sweep_HL(tmin, tmax,s, (tn) 
int tmin, tmax,s; 

int *ltn; 

¢ 

char “sp; 

int 1, lsum, sum,n=2; 

double dt, tn{maxNsweep_change] ; 
tn(0)=ltn£1)=tmin; 
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dot 
if ((++n)>=maxNsweep_change) ¢ 
PutFmt¢ 
"\nToo few H->L sweep changing points. Increase maxNsweep_change=%i", 
maxNsweep_change); 
PutStr("\nPress RETURN to stop and exit"); sp=GetStr(); 
error_flag=SET; 
return(UIO_ERROR); 
> 
(sum=sum; 
sum=(tn(1]; 
dt=( (double)(1.0/tmax)-(double)(1.0/tmin))/(n-1); 
for(i=t;i<n;1++)¢ 
ltn(iti]=nint(tnli)); 
tn(ij=1.0/¢(1.0/tn(i-1)]+dt); 
sum+=nint(tn{i]); 
) 
) while (sum<s); 
itn{--n]+=(s- lsum); 
(tn(0) =n; 
return(UIO_ OK); 
> 


double SweptSquareWaveHceLC ) 
/* Generating alternating swept Square wave starting with high frequency 
down to the low frequency*/ 


if((sweep_count++)==sweep_change_times (Nsweep_change+1] )¢ 
Nsweep_change=(++Nsweep_change)%sweep_change_times (0); 
sweep_count=1; /* renew counting for the next step */ 
one=-one; /* switching state +/- 1 */ 


return(one); 

double SweptSquareWaveL2H(cycle) 

int cycle; 

/* Generating alternating swept square wave starting with low frequency 


up to the high frequency */ 


/* Shaping the spectrum of the swept wave */ 


f f log+=df log; /* Incrementing logarithmic frequency */ 
flog=expl0(fflog); /* Logarithmic frequency */ 
flint=df lin; /* Incrementing linear frequency */ 


df=pow( (double) (mRC- 1-cycle%mRC)/ (double) (mRC-1),3.0)*60.0; 
f=(df*flog+flin)/(df+1.0); 

T+=ts; /* Incrementing Time */ 
return(copysign(1.0,sin(2*PI*f*T))); 

> 


double SweptSineWaveL2H(cycle) 

int cycle; 

/* Generating alternating swept sine wave starting with low frequency 
up to the high frequency */ 


{ 

/* Shaping the spectrum of the swept wave */ 

ff log+=dflos; /* Incrementing logarithmic frequency */ 
flog=expl0(fflog); /* Logarithmic frequency */ 

flint=df lin; /* Incrementing linear frequency */ 


df=pow( (double) (mRC-1-cycleZmRC)/ (double) (mRC-1),3.0)*60.0; 
f=(df*flog+flin)/(df+1.0); 

T+=ts; /* Incrementing Time */ 
return(sin(2*PI*f*T)); 


247 


/* file: window dat.h 


* date: 14 Oct 1991 
ud revised: 14 Oct 1991 
* programmer; Shahar Oror 
* 
* This file contains the definitions for windowing and plotting relative 
* errors of the output layer. 
* 
WHKRKRHHURKKERRRERERRRRRRERERRERRRREEREERRREREEEURER ER EEREREREREREREERERERERREEE | 
#define Ytop 800 /* Topmost window place */ 
#define Ybot 450 /* Bottommost window place */ 
#define Yspace 10 /* Horizontal window spacing */ 
#define Xstart 10 /* Left-most window place */ 
#define Xmsg Xstart+100 /* Left-most message window place */ 
#define Ymsg Ytop /* Bottommost message window place */ 
#define Hmsg 20 /* Height of message window */ 
#define Lmsg 500 /* Length of message window */ 


#define window_length 400 /* Window length */ 

#define window height ((Ytop- Ybot )-(No_of_windows-1)*Yspace)/No_of_windows 
#define max_Nwindows 8  /* maximum 8 of windows */ 

#define msg_ “win 8 /* Message window # */ 

#define max_win_height 75 /* maximum window height */ 

#define Xend (Xstart+winmdow_length) 


#define back_color 0 /* black */ 

#define msg_color 3 /* orange */ 

int No_of_windows=output_lay_ length; /* initial # of windows */ 
int winf lag=RESET; /® flag (SET) indicating windows are open */ 
int Point1Flag{max_Nwindows] =(0); /* first plotting point flag (SET) */ 
double Yscale[{max_Nwindows] ; /* Y-coordinate scale factor */ 
double Yoffset (max_ Nwindows] ; /* Y-coordinate offset */ 


int Line_Color {max _Nwindows] ; 


void open_windows() 

{ 

int wi; /* window index */ 

char *sp,buf (90); 

double min_val,max_val;/*Window limits values*/ 

int xll,yll,/*lower left coordinate*/ 
xur,yur;/*upper right coordinate*/ 


for(wi=0; wi<No_ of windows; wit+){ 
yur=Y top-wi* (window _height+Yspace); 
xur=Xstart+window _length+1; 
xll=Xstart; 
yl l=yur- window | height; 
ug_window(wit1, back color Atl,yll,xur,yur); 
ug_box(wit1,7, 0 ,xlljyll,xur, yur,0); 
sprintf (buf, “\nWindow %1 min_val (0) " (witl)); 
PutStr(buf );sp=GetStr(); 

if(!*sp) sprintf(sp,"0"); 
sscanf(sp,"%f",&min val); 
ug puts(wit1,7,0,1,1,sp,0); 
sprintf(buf,"\nWindow %1 max_val [1] ",(wit1)); 
PutStr(buf);sp=GetStr(); 

1fCl*sp) Sprinti(sp, “1p: 
sscanf(sp,"%f",&max_val); 
ug_puts(witl, 7,0, 1, (window _height-15),sp,0); 
Yscale [wil] =window _height/(max_ val-min_val); 
PutStr(buf); 
Yoffset (wi]=min_val; 
Line_Color (wi] =wi+1; 


) 
winF lag=SET; 
} 
void PlotData(wi, Yval) /* plot next data point */ 
double Yval; 
int wi; 
{ 
int y; /*(relative) y-coordinate */ 
Static int Xold{max_Nwindows] =(0); /*previous X plotting point*/ 
Static int Yold{max_Nwindows] ; /*previous Y plotting point*/ 


/* compute where the point should be plotted */ 
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y=(Yval-Yoffset(wi] )*Yscale(wi]; 
if¢y<1) y=1; 
if(y>=window_ height) y=window_height-1; 


/* clear old cursor out */ 
if(Xold(wi]) ug line(wit+t1,back_color,0,Xold{wi] ,1,Xold(wi] ,window height-1,0); 


/* plot new data point */ 

if(Xold{wi]==0 || Point1Flag(wi] ==RESET) 
ug_point(wit+l,Line Color (wi) ,0,Xxold(wi],y); 

else ug _line(wi+1,Line_Color(wi] ,0,Xold{wi]-1,Yold(wi] ,Xold({wi] ,y,0); 


/* save values for posterity */ 
Yold(wi] =y; 
PointiFlag(wi]=SET; /* indicate that first point already plotted */ 


/* compute next curser location */ 
Xold(wi] ++; 
if (Xold{wi]>=window_length) Xold{wi]=0; 


/* draw new curser line */ 
ug_line(wi+1,15,0,Xold(wi],1,Xold{wi] ,window height-1,0); 
> 
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csv3.1 'file format is Control Strategy Version 3.1 
! Feb 92 (emulbpt.nnc) Back Propagation Network 


| Rae RRR RK KK KKK KK KKK KER RERRRREREKERKRERRRREAREKERRRKRARRRRRRRKER 


i* * 
1* Back Propagation Network Control Strategy, with inupt and output™ 
Lis Regression layers. Suitable for N models (N>=1). 

1* * 
as Modified by: Shahar Oror. February 26th 1992. e 


'* 1.This control strategy allows to test in recall mode, while loading* 
1* the desired output via the ‘read’ function. In Learn mode it works* 


!* identically to the original backprop control strategy. * 
'* Here is an example how to utilize it: * 
i* NINT UIO Read(LayN, mPEs, Datap) * 
1% NINT LayN; - 
1% NINT nPEs; = 
3 SREAL *Datap; . 
sd { . 
* if(LayN==Out_lay) return(UIO_Learn_Output(LayN,nPEs,Datap));* 
os else return( UlO_Learn_Input(LayN,nPEs,Datap)); * 
1* } *® 
1* * 
lie Revision date: 2 March 92 > 


1* 2. The control strategy supports a network with built in delay nea 
1* (regressor) form the output nodes and from the input nodes. 


i* In the learning phase, the desired output (in USRIO) and the * 

'* outputs of the regressors are saved (io Irnrslt for layers 1&2). * 

'* Before the next learning presentation of that same model, the values* 

'* are fed into these layers as if they were there for normal BackProp* 
1* operation. 

1* The values are loaded to the output field of the PE. x 

es Revision date: 31 March 92 - 


1* 3. The network has to have an output buffer layer, and the regress’n* 
1* are taken from the buffer layer (not from the output). This is due* 
1* to NeuralWare implementation (feedback from the output layer always™ 
i* takes the desired value and not the output value. 

\* * 


@eees#s 


i* A similar control strategy is used for training controllers. For* 
1* uniformity, the AUX2 counter comparison was put so that the CS will * 
1* in cases where there is no need to extract information from a middle* 
1* layer. The USRIO program controls the value of AUX2 counter. 

1* 4. In testing the CS supplies "correct" values to the output ee 
1* layer, to check if the learning was ok. In the real game there is a* 
1* direct feedback from the output (practically from the output buffer* 
1* layer), but since the training is not always 100% accurate and the* 
1* recurrent nature of the network, it causes the results to wonder. * 
\* ke 
ty Revision date: 8 April 92 i 
1* 5. Im Recall/Test mode, when a new input is tested, the regressors* 
1* (OutReg & ComReg), that contain the regressed data of the previous* 
1* input that was under test, are initialized. The counter WRK1 is used* 
'* as a switch to determine when to initialize the regressors, and is* 
1* controlled from the USRIO program. 
1 kRakKKKKKkKRKRKKKKkKkKKkkKkKkKkKKkKKKkKkkKkKkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


! 07-Sep-90 Modified to work in learn mode with confusion matrices 
! 


IMASK label op-code operands comment 

L_saR_sa trace aux3 ! 0 set trace option to aux3 

tL saR _sa optclr op:bknc ! 1 do not BKp to PEs w/o conns. BUG!sd 
Li_aR_sa cset recall,0 ! 2 recall count 


| Get input (learn and recall) and desired output (learn only) 
! 


L_saR_sa lset in 3 input layer 

L__-R_sa io read 4 get input data (recall) 
L_saR io lrnin 3 get input data (learn) 
L_saR 


I 
1 
1 

lset cur,1 ! command( input) regression layer(rev.2) 
I 
1 
1 


L_saR_ io lrnin get ComReg(model#) last output(!) data 
L_saR_ lset  cur,1 output regression layer  (rev.2) 
L_saR_ 10 Irnin get OutReg(model#) last output(!) data 


! On the start of test of mew command input, initialize ComReg & OutReg. 
! The counter wrk1 serves as the switch and is controlled from the USRIO. 
' (rev. 5) 

L___R_Sa comp wrki,1 lis wrki<1 
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blt 
lset 
io 
lset 
10 


saR_sa d@dnoreg 


10 
10 
lset 
10 
io 
io 


anoreg Ido not initailize ComReg & OutReg 

In, 1 ! command( input) regression layer 

read linitialize ComReg(model#, input#) (test) 
cur, 1 ! output regression layer 

read tinitialize OutReg(model#, input#) (test) 
lset out,-1 ! output buffer layer 

irnin ' get output(model#)(!) data (learn) 
read ! get output(model#) data(recall,rev. 4) 
out ! 6 output layer 

(rnout ! 7 get desired output (learn) 

read ! 8 get desired output (test) (rev. 1) 
reltst ! 8 get desired output (test) 


Start with the first layer for a forward pass through network 


saR sa 


L 

L 

L 

L 

L 

1 

L 
L_saR_ 
L 

L 

L 

L 

1 

1 

! 

1 

L_saR_ 
L__R_sa 


(set 
ar loop 

math 

lset 


Lemp 
blt 


in ' 9 input layer 
math sumjrnoise|tranjoutput;e=0 ! 10 
sum; [noise/tranj output je=O; fire ! 11 


cur,1 ! 12 next layer 
out ! 13 at output layer ? 
ar loop ! 14 loop till done 


Compute final outputs 


__R_sa 


math 
math 


sumjrnoise;tranjoutputje-=w !15 Recall 
sum| |noise| tran; output je-=wie*=efifire ! 16 learn 


Write results to userio 


! 
Li 
L_saR_ 
! 
! 
4 


L_sar 
L_saR 
L_saR 


Learn cycle - back propagate error. 


lset 
10 
lset 
io 
ccmp 
bit 
lset 
10 
anoe f f 
math 
io 
10 


in,1 'set command( input) regressor layer(rev2 
irnorslt 'save output for next learn of that modl 
cur,1 'set output regressor layer (rev.2) 
irnrslt 'save output for next learn of that modl 
aux2,2 lis aux2<2 (rev. 3) 
anoe f F !'go to output layer (no effector layer) 
IN, aUx2 'set effector layer 
lrnrslt Iwrite its values (for control purposes) 
lset out 'set output layer 

ce=e 1 17 Recall 
write ' 18 recall 
iIrnrslt ' 19 learn 


Store unscaled error in 


current error field. 


Cs 
L_saR_ 
' 

l 

y 

" 


L_saR 
L_saR 
L_saR 
L_saR 
L_sar 
LisaRisa 


It 


lset out ! 20 set output layer 
al Loop math ce=ele*=f'!backpilearnjfire ! 21 
lset  cur,-1 ! 22 previous layer 
lcemp in ! 23 at input layer ? 
bgt al Loop ! 24 loop till done 
trace 0Q ! 25 turn off any trace function 


Zi 
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